Extracting Spectra from NIRCam Grism data

Taking IceAge (ERS-1309; PI: McClure) Data as an Example

Fengwu Sun (Steward Observatory, University of Arizona, 2022)
with significant help from Eiichi Egami (UofA), Nor Pirzkal (STScI) and Marcia Rieke (UofA)
(v0 on July 16, 2022)
(v1 on August 10, 2022)
(v2 on October 21, 2022)


Public PDF of APT File (PID 1309): https://www.stsci.edu/jwst/phase2-public/1309.pdf

A Short Overview:

In this notebook, we will go through the data reduction and spectral extraction pipeline of the JWST/NIRCam Wide-Field Slitless Spectroscopy (WFSS) mode. If you need basic introduction to this observation mode, I would recommend you to take a look of the JDox Page on NIRCam/WFSS. In this notebook, I assume that the potential users have basic knowledge of NIRCam and NIRCam grism, e.g.:

  • There are two independent modules of NIRCam (module A and B);
  • Each module has an instantaneous Field of View (FoV) of ~4.4 arcmin$^2$;
  • There are two dispersion directions called R and C;
  • The dispersion is around 1 nm/pix with a spectral resolution of R~1600 for point sources;
  • When it is taking grism spectra at LW (2.4-5.0µm), the SW channel is taking direct imaging data.

Below is the flow chart for grism data analysis:

flowchart

In this image, all of the items enclosed in red boxes are input fits files you download from MAST, which include (1) direct LW images, (2) direct SW images and (3) LW grism images. Starting from *_rate.fits files (STScI pipeline stage-1 products), you can use direct images in LW and SW to produce final mosaics images and source catalogs (enclosed with green box). These images and source catalogs can be very helpful for your target selections when extracting grism spectra. However, what you actually need to input for grism spectra extraction is (i) astrometric catalog that you want to register SW images onto, and (ii) coordinates (RA and DEC) of targets with interest. As long as you can register the astrometry of the LW grism images, using the spectral tracing models, grism dispersion models and flux calibration functions that I produced, you can extract calibrated 2D grism spectra for any sources you want. We will get to the point on how to do astrometric calibration and spectral calibration later in this notebook.

List of Content

  1. Import necessary modules and define functions

  2. Getting the list of RATE.fits files)

  3. Reduce grism data to level-1.5 (e.g., assign wcs, flatfielding, bkgsubtraction))

  4. Characterize the astrometry of direct imaging in SW

  5. Generating source catalog for each grism exposure

  6. 2D Spectral extraction

  7. Examine 2D & 1D spectra in plots

1. Import necessary python modules

In [24]:
''' Astropy modules (I'm using astropy v5.1) '''
from astropy.stats import sigma_clipped_stats, SigmaClip
from astropy.io import fits, ascii
from astropy.visualization.mpl_normalize import ImageNormalize
from astropy.visualization import MinMaxInterval, PercentileInterval, SqrtStretch, ImageNormalize
from astropy import units as u
from astropy import constants as c
from astropy import wcs, table
from astropy.coordinates import SkyCoord
from astropy.table import Table, Column, vstack
from astropy.time import Time
from astropy.modeling.models import custom_model
from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0 = 70, Om0 = 0.3)

''' Photutils modules (I'm using 1.5.0)'''
from photutils import find_peaks, psf, DAOStarFinder, IRAFStarFinder
from photutils import CircularAperture, aperture_photometry, CircularAnnulus, RectangularAperture, RectangularAnnulus
from photutils import Background2D, MedianBackground, SExtractorBackground
from photutils.segmentation import detect_sources, SourceFinder, SourceCatalog
### photutils background esitmator
sigma_clip = SigmaClip(sigma = 3.)
bkg_estimator = MedianBackground()

''' Numpy and Scipy modules '''
import numpy as np
from scipy import interpolate, integrate, optimize, stats, ndimage, signal

'''Matplotlib modules'''
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patches as patches
from matplotlib.lines import Line2D
import matplotlib.animation as animation

%matplotlib inline

'''Some of the following modules are not used, remove if you cannot import'''
## finds all the pathnames matching a specified pattern
import glob   
## OS commands
import os
## subprocess with shell commands
import subprocess
## get time stamp
import time
## regular expression
import re
## download packadges? I guess I didn't use these:
# import urllib
# import requests
# import pkg_resources
## Belows are used by my legacy codes on grism simulation, maybe no need any more
# import h5py
# import yaml
# import copy

'''You might need some data from mirage simulator: 
       https://github.com/spacetelescope/mirage/
'''
# os.environ["MIRAGE_DATA"] = os.environ["MIRAGE_DATA"] + "/mirage_data/"
'''Import CRDS Server to environmental variable'''
os.environ["CDRS_SERVER_URL"]="https://jwst-cdrs.stsci.edu"

## Multiple processing
from multiprocessing import Pool, Lock
## https://github.com/npirzkal/GRISMCONF
import grismconf
## https://github.com/spacetelescope/pysiaf
import pysiaf 
## https://pypi.org/project/asdf/
import asdf

## if there is any animation to be display with Jupyter notebook:
from IPython.display import HTML

'''Filter some warning'''
import warnings
warnings.filterwarnings("ignore", message="'obsfix' made the change")
warnings.filterwarnings("ignore", message="'datfix' made the change")
warnings.filterwarnings("ignore", message = "Card is too long")

'''Default nircam zeropoint with native pixel size in LW (0.0629 arcsec/pixel)'''
nircam_LW_ZP = -2.5 * np.log10((u.MJy / u.sr * (0.0629*u.arcsec)**2 / (3631 * u.Jy)).cgs.value)
In [25]:
'''
Matplotlib Drawing:
'''
plt.rcParams.update({'font.size': 16})  ## larger font size

def get_corner_pos(ax, loc = 1, edge = 0.02):
    '''get (x, y) for the corner of an axes at given quadrant (1, 2, 3, 4)'''
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    xscale, yscale = ax.get_xscale(), ax.get_yscale()
    if loc not in [1, 2, 3, 4]: raise ValueError('The parameter `loc` must be 1 or 2 or 3 or 4.')
    vec_edge = np.array([1-edge, edge])
    vec_x = np.array([[xlim[1], xlim[0], xlim[0], xlim[1]][loc-1], [xlim[0], xlim[1], xlim[1], xlim[0]][loc-1]])
    vec_y = np.array([[ylim[1], ylim[1], ylim[0], ylim[0]][loc-1], [ylim[0], ylim[0], ylim[1], ylim[1]][loc-1]])
    if xscale == 'linear': x_text = np.dot(vec_edge, vec_x)
    else: x_text = 10 ** np.dot(vec_edge, np.log10(vec_x))
    if yscale == 'linear': y_text = np.dot(vec_edge, vec_y)
    else: y_text = 10 ** np.dot(vec_edge, np.log10(vec_y))
    return (x_text, y_text)

def corner_text(ax, s ='', loc = 1, edge = 0.02, **kwargs):
    '''Add text at the corner of an axes'''
    if loc not in [1, 2, 3, 4]: raise ValueError('The parameter `loc` must be 1 or 2 or 3 or 4.')
    ha = ['right', 'left', 'left', 'right']
    va = ['top', 'top', 'bottom', 'bottom']
    return ax.text(*get_corner_pos(ax, loc, edge), s=s, ha = ha[loc-1], va = va[loc-1], **kwargs)

def rotate(axis, angle):
    '''Fundamental rotation matrices.
    Rotate by angle measured in degrees, about axis 1 2 or 3'''
    if axis not in list(range(1, 4)):
        print ('Axis must be in range 1 to 3')
        return
    r = np.zeros((3, 3))
    ax0 = axis-1 #Allow for zero offset numbering
    theta = angle * np.pi / 180. # radians(angle)
    r[ax0, ax0] = 1.0
    ax1 = (ax0+1) % 3
    ax2 = (ax0+2) % 3
    r[ax1, ax1] = np.cos(theta)
    r[ax2, ax2] = np.cos(theta)
    r[ax1, ax2] = - np.sin(theta)
    r[ax2, ax1] = np.sin(theta)
    return r

def attitude(v2, v3, ra, dec, pa):
    '''This will make a rotation matrix which rotates a unit vector representing a v2, v3 position
    to a unit vector representing an RA, Dec pointing with an assigned position angle
    Described in JWST-STScI-001550, SM-12, section 6.1'''
    # v2, v3 in arcsec, ra, dec and position angle in degrees
    v2d = v2/3600.0
    v3d = v3/3600.0

    # Get separate rotation matrices
    mv2 = rotate(3, -v2d)
    mv3 = rotate(2, v3d)
    mra = rotate(3, ra)
    mdec = rotate(2, -dec)
    mpa = rotate(1, -pa)

    # Combine as mra*mdec*mpa*mv3*mv2
    m = np.dot(mv3, mv2)
    m = np.dot(mpa, m)
    m = np.dot(mdec, m)
    m = np.dot(mra, m)

    return m

'''This is the Grism Dispersion function, will be detailed later'''
def fit_disp_order32(data, 
                     a01, a02, a03, a04, a05, a06, 
                     b01, b02, b03, b04, b05, b06, 
                     c01, c02, c03, #c04, c05, c06,
                     d01, #d02, d03, d04, d05, d06,
                    ):
    ## data is an numpy array of the shape (3, N)
    ##     - data[0]:  x_pixel      --> fit with second-degree polynomial
    ##     - data[1]:  y_pixel      --> fit with second-degree polynomial
    ##     - data[2]:  wavelength   --> fit with third-degree polynomial
    xpix, ypix, dx = data[0] - 1024, data[1] - 1024, data[2] - 3.95
    ## return dx = dx(x_pixel, y_pixel, lambda)
    return ((a01 + (a02 * xpix + a03 * ypix) + (a04 * xpix**2 + a05 * xpix * ypix + a06 * ypix**2)
            ) + 
            (b01 + (b02 * xpix + b03 * ypix) + (b04 * xpix**2 + b05 * xpix * ypix + b06 * ypix**2)
            ) * dx +
            (c01 + (c02 * xpix + c03 * ypix) #+ (c04 * xpix**2 + c05 * xpix * ypix + c06 * ypix**2)
            ) * dx**2 + 
            (d01 #+ (d02 * xpix + d03 * ypix) + (d04 * xpix**2 + d05 * xpix * ypix + d06 * ypix**2)
            ) * dx**3
           ) 

func_fit_wave = fit_disp_order32


'''This is the Spectral Tracing function, will be detailed later'''
def fit_disp_order23(data, 
                     a01, a02, a03, a04, a05, a06, a07, a08, a09, a10,
                     b01, b02, b03, b04, b05, b06, b07, b08, b09, b10,
                     c01, c02, c03, c04, c05, c06, c07, c08, c09, c10,
                    ):
    ## data is an numpy array of the shape (3, N)
    ##     - data[0]:  x_pixel  --> fit with second-degree polynomial
    ##     - data[1]:  y_pixel  --> fit with second-degree polynomial
    ##     - data[2]:  dx     --> fit with second-degree polynomial
    xpix, ypix, dx = data[0] - 1024, data[1] - 1024, data[2]
    ## return dy = dy(x_pixel, y_pixel, d_x)
    return ((a01 + (a02 * xpix + a03 * ypix) + (a04 * xpix**2 + a05 * xpix * ypix + a06 * ypix**2)
             + (a07 * xpix**3 + a08 * xpix**2 * ypix + a09 * xpix * ypix**2 + a10 * ypix**3)
            ) + 
            (b01 + (b02 * xpix + b03 * ypix) + (b04 * xpix**2 + b05 * xpix * ypix + b06 * ypix**2)
             + (b07 * xpix**3 + b08 * xpix**2 * ypix + b09 * xpix * ypix**2 + b10 * ypix**3)
            ) * dx +
            (c01 + (c02 * xpix + c03 * ypix) + (c04 * xpix**2 + c05 * xpix * ypix + c06 * ypix**2)
             + (c07 * xpix**3 + c08 * xpix**2 * ypix + c09 * xpix * ypix**2 + c10 * ypix**3)
            ) * dx**2
           ) 


def grism_conf_preparation(x0 = 1024, y0 = 1024, pupil = 'R', 
                           fit_opt_fit = np.zeros(30), w_opt = np.zeros(16)):
    '''
    Prepare grism configuration, dxs, dys, wavelengths based on input (x0, y0) pixel postion
        and filter/pupil/module information.
    -----------------------------------------------
        Parameters
        ----------  
        x0, y0 : float
            Reference position (i.e., in direct image)
        
        pupil: 'R' or 'C'
            pupil of grism ('R' or 'C')
        
        fit_opt_fit: numpy.ndarray, shape: (30,)
            polynomial parameters in the perpendicular direction, used by function `fit_disp_order23` 
        
        w_opt: numpy.ndarray, shape: (16,)
            polynomial parameters in the dispersed direction, used by function `fit_disp_order32` 
        
        Returns
        -------
        
        dxs, dys : `~numpy.ndarray`
            offset of spectral pixel from the direct imaging position
        
        wavs: `~numpy.ndarray`
            Array of wavelengths corresponding to dxs and dys
    '''
    # Load the Grism Configuration file
    # GConf = grismconf.Config(os.environ['MIRAGE_DATA'] + "/nircam/GRISM_NIRCAM/V3/" +
    #                          "NIRCAM_%s_mod%s_%s.conf" % (filter, module, pupil))
    wave_space = np.arange(2.39, 5.15, 0.01)
    disp_space = func_fit_wave(np.vstack((x0 * np.ones_like(wave_space), y0 * np.ones_like(wave_space), wave_space)), *w_opt)
    inverse_wave_disp = interpolate.UnivariateSpline(disp_space[np.argsort(disp_space)], wave_space[np.argsort(disp_space)], s = 0, k = 1)

    if pupil == 'R':
        dxs = np.arange(int(np.min(disp_space)), int(np.max(disp_space))) - x0%1
        # Compute wavelength of each of the pixels
        wavs = inverse_wave_disp(dxs)
        # Compute the dys values for the same pixels
        dys = fit_disp_order23(np.vstack((x0* np.ones_like(dxs), y0 * np.ones_like(dxs), dxs)), *fit_opt_fit)
        # dxs = np.arange(-1800, 1800, 1) - x0%1
        ## Compute the t values corresponding to the exact offsets (with old grism config)
        # ts = GConf.INVDISPX(order = '+1', x0 = x0, y0 = y0, dx = dxs)
        # dys = GConf.DISPY('+1', x0, y0, ts)
        # wavs = GConf.DISPL('+1', x0, y0, ts)
        # tmp_aper = np.max([0.2, 1.5 * tmp_re_maj * np.cos(tmp_pa), 1.5 * tmp_re_min * np.sin(tmp_pa)])
    elif pupil == 'C':
        # Compute the dys values for the same pixels
        dys = np.arange(int(np.min(disp_space)), int(np.max(disp_space))) - y0%1
        # Compute wavelength of each of the pixels
        wavs = inverse_wave_disp(dys)
        dxs = fit_disp_order23(np.vstack((x0* np.ones_like(dys), y0 * np.ones_like(dys), dys)), *fit_opt_fit)
        # dys = np.arange(-1800, 1800, 1) - y0%1
        ## Compute the t values corresponding to the exact offsets (with old grism config)
        # ts = GConf.INVDISPY(order = '+1', x0 = x0, y0 = y0, dy = dys)
        # dxs = GConf.DISPX('+1', x0, y0, ts)
        # wavs = GConf.DISPL('+1', x0, y0, ts)
    return (dxs, dys, wavs)

'''Lienar Function'''
linear = lambda x, k, b: x * k + b
'''gaussian function'''
gauss = lambda x, x0, flux, FWHM : (flux / FWHM / 1.064467) * np.exp(-(x - x0)**2 / (2 * (FWHM/2.354820)**2))
'''gaussian function + continuum'''
gauss_cont_prof = lambda x, x0, flux, FWHM, k, b : (flux / FWHM / 1.064467) * np.exp(-(x - x0)**2 / (2 * (FWHM/2.354820)**2)) + k * (x - 1024) + b

Not all of the source in your input catalog could yield spectra on your 2D grism image. For example, below is the spectral coverage for F356W grism-R in module A and B:

F356W Grism-R Coverage

In this plot, everything is in the detector frame (X and Y; e.g., the detector range is X=0-2047, Y=0-2047 which is indicated by the red box). Based on pre-launch estimate, the regions that can be perfectly picked off by the POM are shown in shallow grey color, and the coronagraph regions are in white - usally the spectral quality will be affected if your source falls in these regions. Only sources in the blue (module-A) and orange (module-B) can yield complete spectra with the current filter setting (F356W). The pre-launch estimate of full spectral coverage region is indicated by the black boxes, so you can see some differences. The yielded spectra alongside with the dispersion directions are indicated by the arrows. If the source falls in the hatched region, it could yield partial spectral coverage. You can get access to all of the spectral coverage plots of all NIRCam LW filters and grisms at http://u.arizona.edu/~fengwusun/data/spectral_cov/.

The following codes are used to check whether the source could enter Pick-Off Mirror (POM) and yield spectrum:

In [26]:
## module A/B LW POM transmission:
##     These files can be downloaded at: https://github.com/npirzkal/GRISM_NIRCAM/tree/master/V4
modA_LW_POM_fits = fits.open(os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_ModA_trans.fits')
modB_LW_POM_fits = fits.open(os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_ModB_trans.fits')
modA_LW_POM_trans = modA_LW_POM_fits[1].data # np.clip(np.int8(modA_LW_POM_fits[1].data), 0, 1)
modB_LW_POM_trans = modB_LW_POM_fits[1].data # np.clip(np.int8(modB_LW_POM_fits[1].data), 0, 1)
xy_start_POM_modA = modA_LW_POM_fits[1].header['NOMXSTRT'], modA_LW_POM_fits[1].header['NOMYSTRT'] 
xy_start_POM_modB = modB_LW_POM_fits[1].header['NOMXSTRT'], modA_LW_POM_fits[1].header['NOMYSTRT'] 

def is_pickoff(x, y, module = 'A'):
    '''
    Input x/y pixel(s) and module, find whether the source will be picked off by NIRCam POM.
    '''
    ### get module and start point
    if module.upper() == 'A': xy_start, POM_trans = xy_start_POM_modA, modA_LW_POM_trans
    elif module.upper() == 'B': xy_start, POM_trans = xy_start_POM_modB, modB_LW_POM_trans
    else: raise ValueError("Only \"A\" or \"B\" is allowed for module.")
    if type(x) == np.ndarray: tmp_x, tmp_y = np.int32(x), np.int32(y)
    elif type(x) == list: tmp_x, tmp_y = np.array(x, dtype = np.int32), np.array(y, dtype = np.int32)
    elif type(x) == int: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32) 
    elif type(x) == float: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32) 
    else: raise TypeError("Only integers/floats or their list/np.array() are allowed.")
    if len(tmp_x) != len(tmp_y): raise ValueError("x, y list should be of the same length.")
    tmp_x = tmp_x + xy_start[0]
    tmp_y = tmp_y + xy_start[1]
    arr_trans = []
    for j in range(len(tmp_x)):
        if (tmp_x[j] < 0)|(tmp_y[j]<0)|(tmp_x[j] > POM_trans.shape[1]-1)|(tmp_y[j] > POM_trans.shape[0]-1):
            arr_trans.append(0) 
        else:
            arr_trans.append(POM_trans[tmp_x[j], tmp_x[j]])
    return np.array(arr_trans, dtype = int)

def is_pickoff_PS(x, y, module = 'A', filter = 'F322W2', pupil = 'C'):
    '''
    Input x/y pixel(s) and module, 
    find whether the source will be picked off by NIRCam POM and yield partial/complete spectra.
    '''
    ### Get module and start point
    ## partial spectral coverage region fits files (`path_PS1) can be downloaded at: 
    ##        http://gxn.as.arizona.edu/~sunfengwu/spec_cov/
    ## now available for both wide and medium LW filter now.
    path_PS = '/groups/egami/mirage_data/mirage_data/nircam/GRISM_NIRCAM/FSun_SpecCov_%s_%s_%s.fits' % (filter, module, pupil)
    path_POM = os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_Mod%s_trans.fits' % module
    if os.path.isfile(path_PS) == False: 
        raise FileNotFoundError("Spectral Coverage file %s not found!" % path_PS)
    # if os.path.isfile(path_POM) == False: 
    #     raise FileNotFoundError("POM Coverage file %s not found!" %path_POM)
    spec_cov_fits = fits.open(path_PS)
    xy_start = (spec_cov_fits[0].header['NOMXSTRT'], spec_cov_fits[0].header['NOMYSTRT'])
    POM_trans = spec_cov_fits[0].data # * np.int16(fits.getdata(path_POM) == 1)
    spec_cov_fits.close()
    if type(x) == np.ndarray: tmp_x, tmp_y = np.int32(x), np.int32(y)
    elif type(x) == list: tmp_x, tmp_y = np.array(x, dtype = np.int32), np.array(y, dtype = np.int32)
    elif type(x) == int: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32) 
    elif type(x) == float: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32) 
    else: raise TypeError("Only integers/floats or their list/np.array() are allowed.")
    if len(tmp_x) != len(tmp_y): raise ValueError("x, y list should be of the same length.")
    tmp_x = tmp_x + xy_start[0]
    tmp_y = tmp_y + xy_start[1]
    arr_trans = []
    for j in range(len(tmp_x)):
        if (tmp_x[j] < 0)|(tmp_y[j]<0)|(tmp_x[j] > POM_trans.shape[1]-1)|(tmp_y[j] > POM_trans.shape[0]-1):
            arr_trans.append(0) 
        ### For now, I don't trust the spectra from Coronagraphy region (y >~ 2250); 
        ### I'm masking all of these frames:
        elif tmp_y[j] > 2250:
            arr_trans.append(0)
        else:
            arr_trans.append(POM_trans[tmp_y[j], tmp_x[j]])
    return np.array(arr_trans, dtype = float)

2. Get list of RATE files (stage-1 products; from MAST)

After downloading your rate files from MAST https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html or reducing *_uncal.fits to *_rate.fits files, you can get a list of the rate files and start analyzing the observational setups (e.g., filter, module, pupil, targets) of these observations:

In [13]:
'''
List of the rate files on your local directory:
    Here I'm only selecting module A/B long wavelength images by wildcards: `jw*nrc[a-b]long_rate`
'''
# list_rate = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/obs24_F322W2_C/jw*nrc[a-b]long_rate.fits'))
list_rate = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/*/jw*nrc[a-b]long_rate.fits'))


list_rate.sort()

'''List of filter, pupil, module, target, exposure time, etc.'''
list_filter = []
list_pupil = []
list_module = []
list_target = []
list_exp = []
for x in list_rate:
    ### read fits header one by one:
    tmp_hd = fits.getheader(x)
    list_filter.append(tmp_hd['FILTER'])
    list_pupil.append(tmp_hd['PUPIL'])
    list_module.append(tmp_hd['MODULE'])
    list_target.append(tmp_hd['TARGPROP'])
    list_exp.append(tmp_hd['EFFEXPTM'])
list_filter = np.array(list_filter)
list_pupil = np.array(list_pupil)
list_module = np.array(list_module)
list_target = np.array(list_target)
list_exp = np.array(list_exp)
print(np.unique(list_filter), '\n', 
      np.unique(list_pupil), '\n', 
      np.unique(list_module), '\n', 
      np.unique(list_target))
['F322W2' 'F410M' 'F430M' 'F444W'] 
 ['CLEAR' 'GRISMC' 'GRISMR'] 
 ['A' 'B'] 
 ['CHAMMS1-C2-FIELD']

Now at above, you can find out what are the exposures taken in the LW grism (look at the filter and pupil, pupil should be GRISMR or GRISMC). You can also check the number of available exposures using the following command:

In [14]:
### check number of grism exposures
np.sum((list_filter == 'F322W2') & (list_pupil == 'GRISMC') & (list_module == 'A'))
Out[14]:
24

3. Reduce Grism exposures to stage level-1.5 (i.e., bkg subtracted, wcs assigned)

In the following steps, we will reduce the grism data to level 1.5. This includes:

  • World Coordinate System (WCS) assginment.

With the rate files you get above, it usually doesn't contain WCS in the sci fits header. You will need to add WCS information partially using the STScI stage-2 pipeline. The WCS information in such images will be use to convert (RA, Dec) coordinates of your sources to (x, y) arrays assuming the pupil wheel is clear.

  • Flat Fielding

Currently, there is no available NIRCam grism flat field data. One good approximation is to use the image flat field data constructed in-flight, which is quite accurate (guestimate accuracy ~3% per pixel).

  • Background Subtraction

Currently, the theoretical sky background for NIRCam grism observations are all pre-launch estimate and they are not accurate enough. Therefore, I suggest users to construct super sky background using all flat-fielded exposure taken with the same filter+pupil+module combination, and subtract the background there. We might still end up with some background residuals, which could be subtracted using 2D background modeling method.

  • 1/f Noise Subtraction

I'm using the code tshirt/ROEBA which can subtract the horizontal 1/f noise stripe for C grism data. However, right now I don't have good idea on how to deal with R grism data where the dispersion direction is the same as the readout direction, so I'm skipping it.

Below is a plot showing the steps of flat-fielding, background subraction and 1/f noise subtraction on IceAge data (F322W2, module B, grism C):

F322W2 module B grism C reduction (stage-2a)

In [27]:
'''
## Stage 2.1 Background subtraction + hot pixel subtraction
##       2.1.1 get list of rate.fits by module and pupil:   
'''
tmp_filter = 'F444W' # np.unique(list_filter)[0]
list_rate_this_band = list_rate[(list_filter == tmp_filter) & (list_pupil != 'CLEAR')]

print(len(list_rate_this_band), 'files,\n', list_rate_this_band[0])
96 files,
 /xdisk/egami/fengwusun/ERSs/IceAge/all_rate/obs25_F444W_R/jw01309025001_02101_00001_nrcalong_rate.fits
In [28]:
'''Python modules to run JWST pipelines & 1/f noise subtractions'''
from jwst.pipeline import Detector1Pipeline
from jwst import assign_wcs, datamodels, tweakreg
from jwst.flatfield import flat_field
## https://github.com/spacetelescope/crds
import crds 
## https://tshirt.readthedocs.io/en/latest/index.html
from tshirt.pipeline.instrument_specific import rowamp_sub

'''Directory of grism exposures for output -> Change it as your output folder'''
calibrated_dir = '/xdisk/egami/fengwusun/ERSs/IceAge/grism/'
if os.path.isdir(calibrated_dir + 'plots/') == False: 
    os.mkdir(calibrated_dir + 'plots/') # used to save analytic plots
    
def background_grism_stage2(arr_grism_rate, tmp_module, tmp_pupil, tmp_filter):
    '''
    Function to creat Median background for each module & Pupil
    '''
    print('\n>> Create Median Background for %s - module %s - grism %s' % (tmp_filter, tmp_module, tmp_pupil))
    tmp_med_bkg = sigma_clipped_stats(arr_grism_rate, axis = -1, sigma_upper = 2.0, maxiters = 10)[1]
    tmp_med_bkg_path = calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % (tmp_filter, tmp_module, tmp_pupil)
    tmp_med_bkg_fits = fits.HDUList(fits.PrimaryHDU(tmp_med_bkg))
    tmp_med_bkg_fits[0].header['object'] = 'BKG_%s_mod%s_%s' % (tmp_filter, tmp_module, tmp_pupil)
    tmp_med_bkg_fits[0].header['filter'] = tmp_filter
    tmp_med_bkg_fits[0].header['module'] = tmp_module
    tmp_med_bkg_fits[0].header['pupil'] = tmp_pupil
    tmp_med_bkg_fits.writeto(tmp_med_bkg_path, overwrite = True)
    print('>> Saved background models to %s' % tmp_med_bkg_path)
    
def assignwcs_grism_stage2(rate_grism_file):
    ''' 
    Function to reduce grism data to level 1.5 (subtract constructed bkg; assign wcs)
    '''
    rate_grism_fits = fits.open(rate_grism_file)
    tmp_grism_hd = rate_grism_fits[0].header
    ### read distortion assignment
    siaf_file = crds.getreferences(tmp_grism_hd, reftypes = ['distortion'], ignore_cache = False)['distortion']
    ### find available background
    tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']
    print('   read %s' % tmp_path_rate, '(%s - %s - %s)' % (tmp_filter, tmp_module, tmp_pupil))
    ### correct WCS and background (need to fool STScI pipeline to assume it as NRC image)
    grism_wcs_step = assign_wcs.assign_wcs_step.AssignWcsStep(override_distortion = siaf_file)
    rate_grism_fits[0].header['EXP_TYPE'] = 'NRC_IMAGE'
    grism_image = datamodels.image.ImageModel(rate_grism_fits)
    '''don\'t subtract background here'''
    # try:
    #     tmp_bkg_img = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % (tmp_filter, tmp_module, tmp_pupil))
    #     grism_image.data = np.float32(grism_image.data - tmp_bkg_img)
    # except FileNotFoundError:
    #     pass
    grism_with_wcs = grism_wcs_step(grism_image)
    ### correct for flat:
    tmp_flat_path = crds.getreferences({'INSTRUME': tmp_grism_hd['INSTRUME'], 'READPATT': tmp_grism_hd['READPATT'], 
                                        'SUBARRAY': tmp_grism_hd['SUBARRAY'], 
                                        'DATE-OBS': tmp_grism_hd['DATE'].split('T')[0],'TIME-OBS': tmp_grism_hd['DATE'].split('T')[1],
                                        'DETECTOR': tmp_grism_hd['DETECTOR'], 'CHANNEL': tmp_grism_hd['CHANNEL'], 
                                        'MODULE':  tmp_grism_hd['MODULE'], 'EXP_TYPE': 'NRC_IMAGE', #'EXP_TYPE': 'NRC_WFSS',
                                        'FILTER': tmp_grism_hd['FILTER'], 'PUPIL': 'CLEAR', #  'PUPIL': 'GRISMR', # =
                                       })['flat']
    ### error of Flat
    # 1 / np.percentile(fits.getdata(tmp_flat_path, 'err')/ fits.getdata(tmp_flat_path, 'sci'), [16, 50, 84])
    flat_field.do_flat_field(grism_with_wcs, datamodels.FlatModel(tmp_flat_path))
    grism_save_path = rate_grism_file.replace('rate.fits', 'rate_lv1.5.fits')
    grism_save_path = calibrated_dir + grism_save_path.split('/')[-1]
    print("Saving file with proper WCS in", grism_save_path)
    grism_with_wcs.save(grism_save_path)
    
def get_crds_dict_from_fits_header(tmp_hd):
    '''get the CRDS directory from the fits header'''
    crds_dict = {}
    crds_dict['INSTRUME'] = tmp_hd['INSTRUME'].upper()
    crds_dict['READPATT'], crds_dict['SUBARRAY'] = tmp_hd['readpatt'].upper(), tmp_hd['SUBARRAY'].upper()
    crds_dict['DATE-OBS'], crds_dict['TIME-OBS'] = tmp_hd['DATE-OBS'], tmp_hd['TIME-OBS'] 
    crds_dict['DETECTOR'] = tmp_hd['DETECTOR']
    if crds_dict['INSTRUME'] == 'NIRCAM':
        if crds_dict['DETECTOR'] in ['NRCALONG', 'NRCBLONG']: crds_dict['CHANNEL'] = 'LONG'
        else: crds_dict['CHANNEL'] = 'SHORT'
    crds_dict['MODULE'] = tmp_hd['module']
    crds_dict['EXP_TYPE'], crds_dict['FILTER'], crds_dict['PUPIL'] = tmp_hd['EXP_TYPE'], tmp_hd['filter'], tmp_hd['pupil']
    return crds_dict
In [29]:
'''
## Stage 2.2.a WCS corrections + Flat Fielding subtraction
'''
n_procs = 8
### Uncomment the following lines to do this step in parallel
# with Pool(np.min([n_procs, len(list_rate_this_band)])) as pool:
#     pool.map(assignwcs_grism_stage2, list_rate_this_band[:])

### after this step, change the suffix of the files to `rate_lv1.5.fits`
all_rate_list = np.array([calibrated_dir + os.path.basename(x).replace('rate.fits', 'rate_lv1.5.fits')
                          for x in list_rate_this_band])
In [ ]:
'''Check image of flat-fielded grism exposure'''

a_c_sci = fits.getdata(all_rate_list[0])

plt.close()
plt.figure(figsize = (10, 10))
plt.imshow(a_c_sci, 
           vmin = np.nanpercentile(a_c_sci, 0.5), vmax = np.nanpercentile(a_c_sci, 99.5), origin = 'lower')
plt.tight_layout()
In [763]:
'''
## Stage 2.2.b Generate Super Sky Background
'''
print("- Background + Hot pixel construction for %s - " % tmp_filter)
count_ac, count_bc, count_ar, count_br = 0, 0, 0, 0

for i, tmp_path_rate in enumerate(all_rate_list):
    tmp_grism_hd = fits.getheader(tmp_path_rate)
    tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']
    print('>>> [%3d] ' % i, 'read %s' % tmp_path_rate, 
          '(%s %s %s)' % (tmp_filter, tmp_module, tmp_pupil))
    tmp_grism_rate = fits.getdata(tmp_path_rate, 'sci')
    if (tmp_module == 'A') & (tmp_pupil == 'GRISMC'): 
        if count_ac == 0:  arr_grism_rate_ac = tmp_grism_rate
        else: arr_grism_rate_ac = np.dstack((arr_grism_rate_ac, tmp_grism_rate))
        count_ac += 1
    if (tmp_module == 'B') & (tmp_pupil == 'GRISMC'): 
        if count_bc == 0:  arr_grism_rate_bc = tmp_grism_rate
        else: arr_grism_rate_bc = np.dstack((arr_grism_rate_bc, tmp_grism_rate))
        count_bc += 1
    if (tmp_module == 'A') & (tmp_pupil == 'GRISMR'): 
        if count_ar == 0:  arr_grism_rate_ar = tmp_grism_rate
        else: arr_grism_rate_ar = np.dstack((arr_grism_rate_ar, tmp_grism_rate))
        count_ar += 1
    if (tmp_module == 'B') & (tmp_pupil == 'GRISMR'): 
        if count_br == 0:  arr_grism_rate_br = tmp_grism_rate
        else: arr_grism_rate_br = np.dstack((arr_grism_rate_br, tmp_grism_rate))
        count_br += 1
print(count_ar, count_ac, count_br, count_bc)

'''
## generating background by module and pupil: 
'''
### if you only have C grism data, use below:
# list_arr_grism_rate = [arr_grism_rate_ac, arr_grism_rate_bc]
# list_arr_grism_module = ['A', 'B']
# list_arr_grism_pupil = ['GRISMC', 'GRISMC']
# list_arr_grism_filter = [tmp_filter] * 2

### if you only have R grism data, use below:
# list_arr_grism_rate = [arr_grism_rate_ar, arr_grism_rate_br]
# list_arr_grism_module = ['A', 'B']
# list_arr_grism_pupil = ['GRISMR', 'GRISMR']
# list_arr_grism_filter = [tmp_filter] * 2

### if you have both R and C grism data, use below:
list_arr_grism_rate = [arr_grism_rate_ar, arr_grism_rate_ac, arr_grism_rate_br, arr_grism_rate_bc]
list_arr_grism_module = ['A', 'A', 'B', 'B']
list_arr_grism_pupil = ['GRISMR', 'GRISMC', 'GRISMR', 'GRISMC']
list_arr_grism_filter = [tmp_filter] * 4

n_procs = 8
### Uncomment the following lines to do this step in parallel
# with Pool(np.min([n_procs, len(all_rate_list)])) as pool:
#     pool.starmap(background_grism_stage2, zip(list_arr_grism_rate, list_arr_grism_module, 
#                                               list_arr_grism_pupil, list_arr_grism_filter))
- Background + Hot pixel construction for F444W - 
>>> [  0]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00001_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [  1]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00001_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [  2]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00002_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [  3]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00002_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [  4]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00003_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [  5]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00003_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [  6]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00004_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [  7]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00004_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [  8]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00005_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [  9]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00005_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 10]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00006_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 11]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00006_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 12]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00007_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 13]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00007_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 14]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00008_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 15]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00008_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 16]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00009_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 17]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00009_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 18]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00010_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 19]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00010_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 20]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00011_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 21]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00011_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 22]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00012_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 23]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025001_02101_00012_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 24]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00001_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 25]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00001_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 26]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00002_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 27]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00002_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 28]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00003_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 29]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00003_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 30]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00004_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 31]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00004_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 32]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00005_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 33]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00005_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 34]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00006_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 35]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00006_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 36]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00007_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 37]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00007_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 38]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00008_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 39]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00008_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 40]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00009_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 41]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00009_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 42]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00010_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 43]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00010_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 44]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00011_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 45]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00011_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 46]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00012_nrcalong_rate_lv1.5.fits (F444W A GRISMR)
>>> [ 47]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00012_nrcblong_rate_lv1.5.fits (F444W B GRISMR)
>>> [ 48]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00001_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 49]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00001_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 50]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00002_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 51]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00002_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 52]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00003_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 53]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00003_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 54]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00004_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 55]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00004_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 56]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00005_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 57]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00005_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 58]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00006_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 59]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00006_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 60]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00007_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 61]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00007_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 62]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00008_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 63]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00008_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 64]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00009_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 65]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00009_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 66]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00010_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 67]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00010_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 68]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00011_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 69]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00011_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 70]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00012_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 71]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02101_00012_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 72]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00001_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 73]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00001_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 74]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00002_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 75]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00002_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 76]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00003_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 77]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00003_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 78]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00004_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 79]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00004_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 80]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00005_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 81]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00005_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 82]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00006_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 83]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00006_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 84]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00007_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 85]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00007_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 86]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00008_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 87]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00008_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 88]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00009_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 89]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00009_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 90]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00010_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 91]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00010_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 92]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00011_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 93]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00011_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
>>> [ 94]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00012_nrcalong_rate_lv1.5.fits (F444W A GRISMC)
>>> [ 95]  read /xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309026001_02106_00012_nrcblong_rate_lv1.5.fits (F444W B GRISMC)
24 24 24 24

>> Create Median Background for F444W - module A - grism GRISMR

>> Create Median Background for F444W - module A - grism GRISMC

>> Create Median Background for F444W - module B - grism GRISMR

>> Create Median Background for F444W - module B - grism GRISMC
>> Saved background models to /xdisk/egami/fengwusun/ERSs/IceAge/grism/median_bkg_F444W_modA_GRISMR.fits
>> Saved background models to /xdisk/egami/fengwusun/ERSs/IceAge/grism/median_bkg_F444W_modB_GRISMR.fits
>> Saved background models to /xdisk/egami/fengwusun/ERSs/IceAge/grism/median_bkg_F444W_modA_GRISMC.fits
>> Saved background models to /xdisk/egami/fengwusun/ERSs/IceAge/grism/median_bkg_F444W_modB_GRISMC.fits
In [ ]:
'''Check the grism image:'''
tmp_i = np.random.randint(len(all_rate_list))
tmp_i = 51
print(all_rate_list[tmp_i])
tmp_grism_img = fits.getdata(all_rate_list[tmp_i])
tmp_grism_hd = fits.getheader(all_rate_list[tmp_i])
tmp_grism_bkg = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % 
                             (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']))

fig, ax = plt.subplots(2, 2, figsize = (15, 15))
ax = ax.flatten()

'''ax-0: flat-fielded grism image'''
ax[0].imshow(tmp_grism_img,  origin = 'lower', 
             vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[0], s = '(1) Flat-Fielded', color = 'w', loc = 2, weight = 'bold', fontsize = 20)

'''ax-1: try to subtract sigma-clipped median sky background'''
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'): 
    pass ## this part of IceAge data is affected by a few very bright stars
else: tmp_grism_img = tmp_grism_img - tmp_grism_bkg
ax[1].imshow(tmp_grism_img,  origin = 'lower', 
             vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[1], s = '(2) Median Sky BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)

if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
    corner_text(ax[1], s = '(Not applied for modB - GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)
    
# mask-1 for background / 1-over-f noise subtraction: from segment map
_, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:700,100:700])
segment_map = detect_sources(tmp_grism_img - tmp_med, tmp_rms * 2.0, npixels = 100)

### further background subtraction?
sigma_clip = SigmaClip(sigma = 2.5)
bkg_estimator = SExtractorBackground()
if (tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('B', 'R'):
    tmp_mask = segment_map.data!=0
else: tmp_mask = np.zeros_like(segment_map.data)
bkg = Background2D(tmp_grism_img, (24, 24), filter_size = (5, 5), mask = tmp_mask,
                   sigma_clip = sigma_clip, bkg_estimator = bkg_estimator)
tmp_grism_img = tmp_grism_img - bkg.background

'''ax-2: try to subtract modeled 2D sky background'''
ax[2].imshow(tmp_grism_img,  origin = 'lower', 
             vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[2], s = '(3) Modeled 2D BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)

### for grism C: subtract 1/f noise
if tmp_grism_hd['pupil'][-1] == 'C':
    _, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:1700,100:1700].flatten()[::10])
    segment_map = detect_sources(tmp_grism_img, tmp_rms * 5.0, npixels = 20)
    rowSub, modelImg = rowamp_sub.do_backsub(tmp_grism_img, 
                                             amplifiers = 4, backgMask = segment_map.data == 0)
    tmp_grism_img = rowSub

'''ax-3: try to subtract 1/f noise'''
ax[3].imshow(tmp_grism_img,  origin = 'lower', 
             vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[3], s = '(4) 1/f Noise Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
if tmp_grism_hd['pupil'][-1] == 'R':
    corner_text(ax[3], s = '(Not applied for GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)

tmp_title = '%s (%s - %s - %s)' % (os.path.basename(all_rate_list[tmp_i]).replace('_lv1.5.fits', ''),
                                   tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'])
fig.text(0.5, 0.985, s = tmp_title, ha = 'center', va = 'top', fontsize = 24,)
plt.subplots_adjust(left = 0.02, right = 0.98, top = 0.95, bottom = 0.01, wspace = 0.06, hspace = 0.04)

for tmp_ax in ax: 
    tmp_ax.set_xticks([]); tmp_ax.set_yticks([])
# plt.savefig(calibrated_dir + 'plots/%s_mod%s_%s_%2d_%s.pdf' % (tmp_grism_hd['filter'], tmp_grism_hd['module'], 
#                                                                tmp_grism_hd['pupil'], tmp_i, 
#                                                                tmp_title.split(' ')[0]),  dpi = 100)

plt.show()
In [ ]:
'''If necessary, check background / segmentation / 1-over-f noise model'''
plt.figure(figsize = (8, 8))
plt.imshow(bkg.background, origin = 'lower', vmax = np.nanpercentile(bkg.background, 99.5))
print(np.nanpercentile(bkg.background, [0, 100]))
plt.imshow(tmp_mask, origin = 'lower')
# plt.imshow(modelImg)
In [22]:
'''
## Stage 2.2.c Background Ground (1/f noise subtraction) - For loop or pool multi-processing
'''

def my_grism_bkg_subtraction(rate_grism_file):
    print('>>> subtract background for %s ' % rate_grism_file)
    tmp_grism_fits = fits.open(rate_grism_file)
    tmp_grism_img = tmp_grism_fits[1].data
    tmp_grism_hd = tmp_grism_fits[0].header
    tmp_grism_bkg = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % 
                                 (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']))
    plt.close()
    fig, ax = plt.subplots(2, 2, figsize = (15, 15))
    ax = ax.flatten()
    '''ax-0: flat-fielded grism image'''
    ax[0].imshow(tmp_grism_img,  origin = 'lower', 
                 vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
    corner_text(ax[0], s = '(1) Flat-Fielded', color = 'w', loc = 2, weight = 'bold', fontsize = 20)

    '''ax-1: subtract sigma-clipped median sky background'''
    if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'): pass
    else: 
        tmp_grism_img = tmp_grism_img - tmp_grism_bkg
        tmp_grism_fits[0].header['HISTORY'] = 'sigma-clipped median sky background subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d",  time.localtime())
    ax[1].imshow(tmp_grism_img,  origin = 'lower', 
                 vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
    corner_text(ax[1], s = '(2) Median Sky BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
    if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
        corner_text(ax[1], s = '(Not applied for modB - GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)

    # mask-1: from segment map
    _, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:700,100:700])
    segment_map = detect_sources(tmp_grism_img - tmp_med, tmp_rms * 1.0, npixels = 100)

    ### further background subtraction?
    sigma_clip = SigmaClip(sigma = 2.5)
    bkg_estimator = SExtractorBackground()
    if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
        tmp_mask = segment_map.data!=0
#         bkg = Background2D(tmp_grism_img, (36, 36), filter_size = (5, 5), mask = tmp_mask,
#                            sigma_clip = sigma_clip, bkg_estimator = bkg_estimator, exclude_percentile=50.)
    else: 
        tmp_mask = np.zeros_like(segment_map.data)

    bkg = Background2D(tmp_grism_img, (24, 24), filter_size = (5, 5), mask = tmp_mask,
                       sigma_clip = sigma_clip, bkg_estimator = bkg_estimator, exclude_percentile=50.)
    tmp_grism_img = tmp_grism_img - bkg.background
    tmp_grism_fits[0].header['HISTORY'] = 'SExtractor-modeled 2D background subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d",  time.localtime())
    '''ax-2: subtract modeled 2D sky background'''
    ax[2].imshow(tmp_grism_img,  origin = 'lower', 
                 vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
    corner_text(ax[2], s = '(3) Modeled 2D BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)

    ### for grism C: subtract 1/f noise
    if tmp_grism_hd['pupil'][-1] == 'C':
        _, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:1700,100:1700].flatten()[::10])
        segment_map = detect_sources(tmp_grism_img, tmp_rms * 5.0, npixels = 20)
        rowSub, modelImg = rowamp_sub.do_backsub(tmp_grism_img, 
                                                 amplifiers = 4, backgMask = segment_map.data == 0)
        tmp_grism_img = rowSub
        tmp_grism_fits[0].header['HISTORY'] = '1/f noise subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d",  time.localtime())

    '''ax-3: subtract 1/f noise'''
    ax[3].imshow(tmp_grism_img,  origin = 'lower', 
                 vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
    corner_text(ax[3], s = '(4) 1/f Noise Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
    if tmp_grism_hd['pupil'][-1] == 'R':
        corner_text(ax[3], s = '(Not applied for GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)

    tmp_title = '%s (%s - %s - %s)' % (os.path.basename(rate_grism_file).replace('_lv1.5.fits', ''),
                                       tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'])
    fig.text(0.5, 0.985, s = tmp_title, ha = 'center', va = 'top', fontsize = 24,)
    plt.subplots_adjust(left = 0.02, right = 0.98, top = 0.95, bottom = 0.01, wspace = 0.06, hspace = 0.04)

    for tmp_ax in ax: 
        tmp_ax.set_xticks([]); tmp_ax.set_yticks([])
    plt.savefig(calibrated_dir + 'plots/%s_mod%s_%s_%s.pdf' % (tmp_grism_hd['filter'], tmp_grism_hd['module'], 
                                                               tmp_grism_hd['pupil'], tmp_title.split(' ')[0]),  dpi = 100)
    plt.close()
    
    '''Write Fits File'''
    tmp_grism_fits[1].data = tmp_grism_img
    tmp_grism_fits.writeto(rate_grism_file, overwrite = True)
    tmp_grism_fits.close()

n_procs = 8

### Uncomment the following lines to do this step in parallel:
# with Pool(np.min([n_procs, len(all_rate_list)])) as pool:
#     pool.map(my_grism_bkg_subtraction, all_rate_list)

### Or, uncomment the following lines to do this step in loop:
# for i, x in enumerate(all_rate_list):
#     print(i); my_grism_bkg_subtraction(x)
    
    

At this step, you should have completed the stage-2a reduction for the LW grism images. The files (jw*nrc[a-b]long_rate_lv1.5.fits) you get are flat-fielded, wcs-assigned and background subtracted LW grism images, almost ready for spectral extraction. However, there is one more thing that you will need to pay attention, which is astrometry.

The spectral and spatial resolution of NIRCam (grism) is so high that a small astrometric offset (0.3" -> 5 LW pix) can give you a completely trash 2D or 1D spectra. The quality of astrometric correction usually limits the quality of spectral extraction and flux/wavelength calibration. One good news is that we have SW simultaneous direct imaging data which can be used to empirically calibrate the LW grism astrometry. This is, assuming that when SW direct imaging data is registered to certain astrometric catalog by applying dRA, dDEc and theta (i.e., x/y offset and rotation), the same offset will be applied to the LW grism fits header and the LW grism astrometry should be good. Such a registration relies on the internal alignment between NIRCam SW and LW instrument aperture and it may introduce astrometric residuals. However, these residuals should be stable and have been included in the grism spectral tracing models.

When using the different SW/LW filter combinations, different filters will also introduce astrometric offset called "filter offset". In my personal opinion, this has not been well understood and calibrated throughout commissioning and it can cause most of the tracing and wavelength error (a few pixels). We will get to this point later.

4. Astrometry Characterization using SW Direct imaging

In [32]:
'''Getting all SW rate fits files by wildcards'''
list_rate_sw = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/*/jw*_nrc[a-b][1-4]_rate.fits'))

list_rate_sw.sort()
print(len(list_rate_sw))

'''List of filter, pupil, module, target, exposure time, etc.'''

list_sw_filter = []
list_sw_pupil = []
list_sw_module = []
list_sw_target = []
list_sw_exp = []
for x in list_rate_sw:
    tmp_hd = fits.getheader(x)
    list_sw_filter.append(tmp_hd['FILTER'])
    list_sw_pupil.append(tmp_hd['PUPIL'])
    list_sw_module.append(tmp_hd['MODULE'])
    list_sw_target.append(tmp_hd['TARGPROP'])
    list_sw_exp.append(tmp_hd['EFFEXPTM'])
list_sw_filter = np.array(list_sw_filter)
list_sw_pupil = np.array(list_sw_pupil)
list_sw_module = np.array(list_sw_module)
list_sw_target = np.array(list_sw_target)

### print SW filters and modules
print(np.unique(list_sw_filter), np.unique(list_sw_module), np.unique(list_sw_target))
944
['F140M' 'F150W' 'F182M' 'F200W'] ['A' 'B'] ['CHAMMS1-C2-FIELD']

For IceAge data in our example, if I'm looking at the F444W grism data, the SW simultaneous exposures are always taken in F200W filter, so I'm only selecting SW rate files with list_sw_filter = 'F200W' . For these data, I have run the following steps in another jupyter notebook:

  • reduced all *_rate.fits files to stage-2 as *_cal_fits, including flat-fielding, distortion correction, photometric zeropoint assignment and 1/f noise subtraction.

  • register the WCS of all *_cal_fits to an astrometric catalog, constructed at LW.

Therefore, it might be good to detect stars in stage-2 products *_cal_fits, but the WCS information there has been corrected. What we will need is catalogs of stars in all SW images without any astrometric correction, so we can correct for x/y offset and rotration ourselves and pass that to LW. If you haven't corrected anything for _rate.fits and _cal.fits files, great, you can directly run daofind() on the _cal.fits files, get catalogs and compute the astrometric offsets. If you have corrected the astrometry (like me, I did SW/LW imaging first for all IceAge data), we will restore the original astrometry to extract the catalog. Therefore, for the following code block, please modify the scripts to let it suit your cases.

In [34]:
all_rate_sw_list = list_rate_sw[list_sw_filter == 'F200W']
print(len(all_rate_sw_list), 'SW rate files selected.')

'''Path to my direct imaging folder 
 which has rate files, wcs uncorrected (original),
        and cal files, wcs corrected   (modified).
'''
direct_image_dir = '/xdisk/egami/fengwusun/ERSs/IceAge/direct/'

# i = 72
# tmp_rate_sw = all_rate_sw_list[i]

def my_daofind_sw_fits(tmp_rate_sw):
    '''
    Since SW images have been astrometrically corrected and passed through stage-2 pipelines,
    Here we only restore the original atrometry and use the images provided by stage-2 pipeline
    '''
    ## assign wcs 
    rate_sw_fits = fits.open(tmp_rate_sw)
    rate_sw_hd = rate_sw_fits[0].header
    ### read distortion assignment
    siaf_file = crds.getreferences(rate_sw_hd, reftypes = ['distortion'], ignore_cache = False)['distortion']
    ### assign default WCS
    rate_sw_wcs_step = assign_wcs.assign_wcs_step.AssignWcsStep(override_distortion = siaf_file)
    rate_sw_IM = datamodels.image.ImageModel(rate_sw_fits)
    rate_sw_IM_wcs = rate_sw_wcs_step(rate_sw_IM)
    tmp_wcs = rate_sw_IM_wcs.get_fits_wcs()       ## this is the default wcs from observational parameters
    rate_sw_IM_wcs.close()
    rate_sw_fits.close()

    ### well-processed SW cal.fits files 
    direct_sw_fits = fits.open(direct_image_dir + os.path.basename(tmp_rate_sw).replace('_rate.fits', '_cal.fits'))
    tmp_coord_center = wcs.utils.pixel_to_skycoord(1024.5, 1024.5, tmp_wcs)
    ### construct detection image
    tmp_detect_img = np.nan_to_num(direct_sw_fits['sci'].data)
    bkg_estimator = MedianBackground()
    bkg = Background2D(tmp_detect_img, (64, 64), filter_size = (5, 5),
                       sigma_clip = SigmaClip(sigma = 3.), bkg_estimator = bkg_estimator)
    tmp_detect_std =  sigma_clipped_stats(tmp_detect_img[100:1500,100:1500].flatten()[::7], 
                                          sigma = 3, maxiters = 10)[-1]
    tmp_detect_img = (tmp_detect_img - bkg.background) / tmp_detect_std
    # tmp_detect_img[direct_sw_fits['dq'].data >= 1] = 0
    
    ### run DAOStarFinder - change parameters (FWHM or Threshold if needed)
    daofind = DAOStarFinder(fwhm = 2.5, threshold = 8.)
    tmp_tb_daofind = daofind(tmp_detect_img)
    tmp_tb_daofind['detector'] = direct_sw_fits[0].header['detector'].lower()
    tmp_tb_daofind['skycoord'] = wcs.utils.pixel_to_skycoord(tmp_tb_daofind['xcentroid'], tmp_tb_daofind['ycentroid'], tmp_wcs)
    ### Save SW catalogs 
    ascii.write(tmp_tb_daofind, calibrated_dir + 'astrom/%s' % os.path.basename(tmp_rate_sw).replace('_rate.fits', '_daofind.cat'),
                format = 'ecsv', overwrite = True)
    direct_sw_fits.close()

if os.path.isdir(calibrated_dir + 'astrom/') == False: 
    os.mkdir(calibrated_dir + 'astrom/') # used to save astrometry-related files
    
n_procs = 8
### Uncomment the following lines to generate source catalogs in parallel:
# with Pool(np.min([n_procs, len(all_rate_sw_list)])) as pool:
#     pool.map(my_daofind_sw_fits, all_rate_sw_list)

'''If some cores are corrupted, run the following for unfinished fits'''
idx_unfinished = np.array([os.path.isfile(calibrated_dir + 'astrom/%s' % os.path.basename(x).replace('_rate.fits', '_daofind.cat')) == False
                           for x in all_rate_sw_list])
print(np.sum(idx_unfinished), 'unfinished.')
# with Pool(np.min([n_procs, len(all_rate_sw_list[idx_unfinished])])) as pool:
#     pool.map(my_daofind_sw_fits, all_rate_sw_list[idx_unfinished])
384 SW rate files selected.
0 unfinished

With the SW catalogs produced above, you can then use an astrometric catalog to calibrate the astrometry of each SW exposures. This information will be recorded in a table called tb_sw_astrometry and will be passed to the next step for LW grism astrometry correction.

Read LW direct imaging catalog, compare with SW extracted catalog for grism WCS correction

In this example, the astrometric catalog for IceAge data is the LW direct imaging catalog in F430M band. You can get access to it at http://gxn.as.arizona.edu/~sunfengwu/IceAge/mosaics/daofind_IceAge_NIRCam_F430M_source.cat.

In [35]:
'''Read astrometric catalog'''
tb_LW = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F430M_source.cat')
tb_LW = tb_LW[tb_LW['flags'] == 0]
tb_LW[:3]
Out[35]:
Table length=3
idxcentroidycentroidsharpnessroundness1roundness2npixskypeakfluxmagRADECflux_r2flux_r4fluxerr_r4dao_indicatoraper_indicatorSN_r4flagsflux_r4_corrfluxerr_r4_corrabmagabmag_err
int64float64float64float64float64float64int64float64float64float64float64float64float64float64float64float64float64float64float64int64float64float64float64float64
5215901.851452.680.5763-0.1403-0.0955250.0847.323142.092-5.381166.2459073-77.42661715144.65537281.72940.4157-0.01360.04073494.03208931.274911.720416.70420.0014
4981458.421428.9130.4475-0.08040.1706250.0589.91796.235-4.958166.5861082-77.42705484310.22016295.2860.3964-0.04040.07393167.64907721.370511.383316.86220.0016
8054144.0711986.8860.5002-0.2591-0.0293250.0636.293106.253-5.066166.3805331-77.41778154230.28626244.43850.1729-0.01770.08547204.58107659.00444.384616.8710.0006
In [37]:
'''*_cal.fits files in this band (F200W)'''
list_img_cal_this_band = list_rate_sw[list_sw_filter == 'F200W']
list_img_cal_this_band = np.array([calibrated_dir + '../direct/' + os.path.basename(x).replace('_rate.fits', '_cal.fits') 
                                   for x in list_img_cal_this_band])

In the following code block, you will run a loop to compute the x/y offset and rotation of four SW exposures associated with the LW grism exposure (e.g., A1-4 for module A) and save the astrometric shift information in the tb_sw_astrometry table. Note that module A/B might have different offsets so you have to group them separately.

In [38]:
'''SW astrometry table'''
tb_sw_astrometry = Table(names = ['expName','N_match', 'dRA', 'dDEC', 'dRA_err', 'dDEC_err', 'theta', 'theta_err'], 
                         dtype = ['S40', 'i4', 'f4', 'f4' , 'f4', 'f4', 'f4', 'f4'])
tb_sw_astrometry.meta['comments'] = ['RA/Dec offsets are in arcsec and absolute', 
                                     'i.e., cos(Dec) projection has been considered.',
                                     '  >> dRA  =  RA_daofind -  RA_gaia',
                                     '  >> dDEC = DEC_daofind - DEC_gaia',
                                     'generated by F. Sun (steward) on %s' % time.strftime("%Y/%m/%d",  time.localtime())]
for colname in tb_sw_astrometry.colnames[2:]: tb_sw_astrometry[colname].info.format = '.4f'
sigma_clip = SigmaClip(sigma = 2.)

'''Run the following loop to compute astrometric offsets for each group of SW exposures.'''
print('>>>  Astrometry (RA/DEC) offsets calculation for %3d SW Frames' % len(list_img_cal_this_band))
for k in range(len(list_img_cal_this_band) // 4):
    ## For each cal.fits, get all the paths of csal.fits from all four detector.
    tmp_img_cal_path = list_img_cal_this_band[k * 4]
    tmp_img_cal_path_base = os.path.basename(tmp_img_cal_path)[:30]
    tmp_obsnum = int(tmp_img_cal_path_base[7:10])
    tmp_grism_LW_path = calibrated_dir + tmp_img_cal_path_base + 'long_rate_lv1.5.fits'
    tmp_grism_hd_sci = fits.getheader(tmp_grism_LW_path, 1)
#     if tmp_obsnum not in [125, 126]: continue
    print('%3d >> %s ' % (k, tmp_img_cal_path_base))
    tmp_img_cal_hd = fits.getheader(tmp_img_cal_path, 0)
    tmp_detector = tmp_img_cal_hd['DETECTOR'].lower()
    tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
    tmp_time_obs = Time(tmp_img_cal_sci_hd['MJD-AVG'], format ='mjd')
    tmp_detectors = np.array(['%s%d' % (tmp_detector[:-1], tmp_det) for tmp_det in [1, 2, 3, 4]])
    tmp_img_cal_paths = np.array([tmp_img_cal_path.replace(tmp_detector, x) for x in tmp_detectors])
    
    
    
    ## Run DAOFIND to get a star catalog:
    daofind = DAOStarFinder(fwhm = 3.0, threshold = 8.)
    for j, tmp_img_cal_path in enumerate(tmp_img_cal_paths):
        tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
        '''Directly use prepared SW DAOFIND catalog'''
        tmp_tb_daofind = ascii.read(calibrated_dir + 'astrom/%s' % os.path.basename(tmp_img_cal_path).replace('_cal.fits', '_daofind.cat').replace('_rate.fits', '_daofind.cat'),
                                    format = 'ecsv')
        '''no need to do daofind again'''
#         tmp_wcs = wcs.WCS(tmp_img_cal_sci_hd)
#         tmp_coord_center = SkyCoord(tmp_img_cal_sci_hd['crval1'], tmp_img_cal_sci_hd['crval2'], unit = (u.deg, u.deg))
#         tmp_detect_img = np.nan_to_num(fits.getdata(tmp_img_cal_path, 'sci'))
#         bkg_estimator = MedianBackground()
#         bkg = Background2D(tmp_detect_img, (64, 64), filter_size = (5, 5),
#                            sigma_clip = SigmaClip(sigma = 3.), bkg_estimator = bkg_estimator)
#         tmp_detect_std =  sigma_clipped_stats(tmp_detect_img[100:500,100:500], sigma = 3, maxiters = 10)[-1]
#         tmp_detect_img = (tmp_detect_img - bkg.background) / tmp_detect_std
#         # tmp_detect_img[fits.getdata(tmp_img_cal_path, 'dq') >= 1] = 0
#         # tmp_detect_img = np.nan_to_num(fits.getdata(tmp_img_cal_path, 'sci')/np.clip(fits.getdata(tmp_img_cal_path, 'err'), 1e-8, 1e8))
#         tmp_tb_daofind = daofind(tmp_detect_img)
#         tmp_tb_daofind['detector'] = tmp_detectors[j]
#         tmp_tb_daofind['skycoord'] = wcs.utils.pixel_to_skycoord(tmp_tb_daofind['xcentroid'], tmp_tb_daofind['ycentroid'], tmp_wcs)
        if j == 0 : tb_daofind = tmp_tb_daofind
        else: tb_daofind = vstack((tb_daofind, tmp_tb_daofind))
    tb_daofind = tb_daofind[np.argsort(tb_daofind['peak'])][-200:]
    tmp_coord_daofind = tb_daofind['skycoord'] 
    
    
    '''Astrometric Source RA & DECs'''
    # tmp_coord_gaia = SkyCoord(tmp_RA, tmp_DEC, unit = (u.deg, u.deg))
    # arg_LW_bright = np.where((tb_LW['MAG_AUTO'].data < 22.5) & (tb_LW['MAG_AUTO'].data > 16.5))[0]
    # tmp_coord_LW = SkyCoord(tb_LW['X_WORLD'][arg_LW_bright], tb_LW['Y_WORLD'][arg_LW_bright], 
    #                           unit = (u.deg, u.deg))
    arg_LW_bright = np.where((tb_LW['abmag'].data < 22.5) & (tb_LW['abmag'].data > 16.5))[0]
    tmp_coord_LW = SkyCoord(tb_LW['RA'][arg_LW_bright], tb_LW['DEC'][arg_LW_bright], 
                              unit = (u.deg, u.deg))
    ### rough center of the frame
    tmp_coord_center = SkyCoord(np.median(tmp_coord_daofind.ra), np.median(tmp_coord_daofind.dec))
    ### only select sources close to the center
    tmp_coord_LW = tmp_coord_LW[tmp_coord_center.separation(tmp_coord_LW) < 4 * u.arcmin]
    
    ## Cross match DAOFind Catalog with Astrometric Catalog
    idx_daofind, d2d, _ = tmp_coord_LW.match_to_catalog_sky(tmp_coord_daofind)
    idx_LW = np.where(d2d < 2.5 * u.arcsec)[0]  ### Change the criterion if necessary
    idx_daofind = idx_daofind[idx_LW]

    ## Compute RA/DEC offset
    tmp_ra_offset = ((tmp_coord_daofind[idx_daofind].ra - tmp_coord_LW[idx_LW].ra) 
                     * np.cos(tmp_coord_LW[idx_LW].dec)).to(u.arcsec).value
    tmp_dec_offset = (tmp_coord_daofind[idx_daofind].dec - tmp_coord_LW[idx_LW].dec).to(u.arcsec).value

    ## sigma-clipped RA/DEC offset
    tmp_ra_offset_clipped = sigma_clip(tmp_ra_offset)
    tmp_dec_offset_clipped = sigma_clip(tmp_dec_offset)
    arg_clipped = np.where((tmp_ra_offset_clipped.mask == False) & (tmp_dec_offset_clipped.mask == False))[0]
    tmp_ra_offset_clipped, tmp_dec_offset_clipped = tmp_ra_offset[arg_clipped], tmp_dec_offset[arg_clipped]

    tmp_ra_offset_med, tmp_ra_offset_std = np.median(tmp_ra_offset_clipped), np.std(tmp_ra_offset_clipped)
    tmp_dec_offset_med, tmp_dec_offset_std = np.median(tmp_dec_offset_clipped), np.std(tmp_dec_offset_clipped)
    
    ## Compute rotation:
    poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec, ydata = tmp_ra_offset)[0]
    poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra, ydata = tmp_dec_offset)[0]
    ## iter-1
    arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
    poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])[0]
    arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
    poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])[0]
    ## iter-2
    arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
    poly_dec_dRA, pcov_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])
    arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
    poly_ra_dDEC, pccov_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])
    perr_dec_dRA = np.sqrt(np.diag(pcov_dec_dRA))
    perr_ra_dDEC = np.sqrt(np.diag(pccov_ra_dDEC))
    ## combine two direction
    cos_factor = np.cos(np.deg2rad(np.median(tmp_coord_daofind[idx_daofind].dec.value)))
    theta_as_deg = (poly_dec_dRA[0] / perr_dec_dRA[0] - poly_ra_dDEC[0] / perr_ra_dDEC[0]) / (1/perr_dec_dRA[0] + cos_factor / perr_ra_dDEC[0]) / 3600.
    theta_as_deg_err = np.sum(perr_dec_dRA[0]**2 + (perr_ra_dDEC[0] / cos_factor)**2)**0.5 / 2. / 3600.
    theta_as_deg = np.rad2deg(theta_as_deg)
    theta_as_deg_err = np.rad2deg(theta_as_deg_err)
    
    tmp_ra_offset_med = tmp_grism_hd_sci['CRVAL2'] * poly_dec_dRA[0] + poly_dec_dRA[1]
    tmp_dec_offset_med = tmp_grism_hd_sci['CRVAL1'] * poly_ra_dDEC[0] + poly_ra_dDEC[1]

    tmp_ra_offset_std = sigma_clipped_stats(tmp_ra_offset)[2]
    tmp_dec_offset_std = sigma_clipped_stats(tmp_dec_offset)[2]
    print(' RA_image -  RA_ref = %.3f"±%.3f"' % (tmp_ra_offset_med,  tmp_ra_offset_std  / np.sqrt(len(tmp_ra_offset)) ))
    print('DEC_image - DEC_ref = %.3f"±%.3f"' % (tmp_dec_offset_med, tmp_dec_offset_std / np.sqrt(len(tmp_ra_offset)) ))
    print('     rotation theta = %.3f±%.3f deg' % (theta_as_deg, theta_as_deg_err))
    tb_sw_astrometry.add_row([tmp_img_cal_path_base, len(arg_clipped),
                              tmp_ra_offset_med, tmp_dec_offset_med, 
                              tmp_ra_offset_std, tmp_dec_offset_std,
                              theta_as_deg, theta_as_deg_err])
    plt.close()
    plt.subplots(1, 1, figsize = (5, 5))
    plt.plot(tmp_ra_offset, tmp_dec_offset, marker = '.', color = 'gray', ls = 'none')
    plt.plot(tmp_ra_offset_clipped, tmp_dec_offset_clipped, marker = '.', color = 'k', ls = 'none')
    plt.plot(tmp_ra_offset_med, tmp_dec_offset_med, marker = '+', color = 'red', ms = 20, mew = 2, ls = 'none')
    plt.axhline(0, color = 'grey', ls = '--')
    plt.axvline(0, color = 'grey', ls = '--')
    plt.gca().set(xlabel = r'$\Delta$RA', ylabel = r'$\Delta$Dec')
    if np.max(np.abs(np.array([tmp_ra_offset_med, tmp_dec_offset_med]))) > 1:
        plt.gca().set(aspect = 1, xlim = (-1.5, 1.5), ylim = (-1.5, 1.5), 
                      xticks = np.arange(-1.5, 1.51, 0.5), yticks = np.arange(-1.5, 1.51, 0.3))
    elif  np.max(np.abs(np.array([tmp_ra_offset_med, tmp_dec_offset_med]))) > 0.5:
        plt.gca().set(aspect = 1, xlim = (-1., 1.), ylim = (-1., 1.), 
                      xticks = np.arange(-0.8, 1.1, 0.4), yticks = np.arange(-1., 1.1, 0.2))
    else:
        plt.gca().set(aspect = 1, xlim = (-0.5, 0.5), ylim = (-0.5, 0.5), 
                      xticks = np.arange(-0.4, 0.51, 0.2), yticks = np.arange(-0.5, 0.51, 0.1))
                  
    corner_text(plt.gca(), s = ' RA_image -  RA_ref = %.3f"±%.3f"' % (tmp_ra_offset_med,  tmp_ra_offset_std / np.sqrt(len(tmp_ra_offset)) ),
                loc = 1, fontsize = 12)
    corner_text(plt.gca(), s = 'DEC_image - DEC_ref = %.3f"±%.3f"' % (tmp_dec_offset_med, tmp_dec_offset_std/ np.sqrt(len(tmp_ra_offset)) ),
                loc = 4, fontsize = 12)
    plt.title(tmp_img_cal_path_base, fontsize = 14, )
    plt.tight_layout()
    # plt.savefig('/home/u24/fengwusun/IceAge/SW_astrometry/%s_astrom.jpg' % tmp_img_cal_path_base, dpi = 100)
    # break
>>>  Astrometry (RA/DEC) offsets calculation for 384 SW Frames
  0 >> jw01309025001_02101_00001_nrca 
 RA_image -  RA_ref = -0.387"±0.009"
DEC_image - DEC_ref = 0.210"±0.009"
     rotation theta = 0.100±0.003 deg
  1 >> jw01309025001_02101_00001_nrcb 
 RA_image -  RA_ref = -0.276"±0.013"
DEC_image - DEC_ref = -0.057"±0.008"
     rotation theta = 0.095±0.003 deg
  2 >> jw01309025001_02101_00002_nrca 
 RA_image -  RA_ref = -0.387"±0.009"
DEC_image - DEC_ref = 0.211"±0.008"
     rotation theta = 0.101±0.003 deg
  3 >> jw01309025001_02101_00002_nrcb 
 RA_image -  RA_ref = -0.272"±0.014"
DEC_image - DEC_ref = -0.056"±0.008"
     rotation theta = 0.094±0.003 deg
  4 >> jw01309025001_02101_00003_nrca 
 RA_image -  RA_ref = -0.386"±0.009"
DEC_image - DEC_ref = 0.208"±0.009"
     rotation theta = 0.099±0.003 deg
  5 >> jw01309025001_02101_00003_nrcb 
 RA_image -  RA_ref = -0.270"±0.013"
DEC_image - DEC_ref = -0.060"±0.008"
     rotation theta = 0.096±0.003 deg
  6 >> jw01309025001_02101_00004_nrca 
 RA_image -  RA_ref = -0.386"±0.009"
DEC_image - DEC_ref = 0.209"±0.009"
     rotation theta = 0.098±0.003 deg
  7 >> jw01309025001_02101_00004_nrcb 
 RA_image -  RA_ref = -0.273"±0.013"
DEC_image - DEC_ref = -0.059"±0.008"
     rotation theta = 0.095±0.003 deg
  8 >> jw01309025001_02101_00005_nrca 
 RA_image -  RA_ref = -0.400"±0.010"
DEC_image - DEC_ref = 0.215"±0.008"
     rotation theta = 0.098±0.003 deg
  9 >> jw01309025001_02101_00005_nrcb 
 RA_image -  RA_ref = -0.281"±0.012"
DEC_image - DEC_ref = -0.052"±0.008"
     rotation theta = 0.093±0.003 deg
 10 >> jw01309025001_02101_00006_nrca 
 RA_image -  RA_ref = -0.402"±0.009"
DEC_image - DEC_ref = 0.216"±0.009"
     rotation theta = 0.098±0.003 deg
 11 >> jw01309025001_02101_00006_nrcb 
 RA_image -  RA_ref = -0.286"±0.012"
DEC_image - DEC_ref = -0.053"±0.008"
     rotation theta = 0.097±0.003 deg
 12 >> jw01309025001_02101_00007_nrca 
 RA_image -  RA_ref = -0.400"±0.010"
DEC_image - DEC_ref = 0.219"±0.008"
     rotation theta = 0.098±0.003 deg
 13 >> jw01309025001_02101_00007_nrcb 
 RA_image -  RA_ref = -0.283"±0.013"
DEC_image - DEC_ref = -0.049"±0.008"
     rotation theta = 0.094±0.003 deg
 14 >> jw01309025001_02101_00008_nrca 
 RA_image -  RA_ref = -0.400"±0.009"
DEC_image - DEC_ref = 0.214"±0.009"
     rotation theta = 0.098±0.003 deg
 15 >> jw01309025001_02101_00008_nrcb 
 RA_image -  RA_ref = -0.284"±0.012"
DEC_image - DEC_ref = -0.055"±0.008"
     rotation theta = 0.095±0.003 deg
 16 >> jw01309025001_02101_00009_nrca 
 RA_image -  RA_ref = -0.418"±0.012"
DEC_image - DEC_ref = 0.225"±0.009"
     rotation theta = 0.098±0.003 deg
 17 >> jw01309025001_02101_00009_nrcb 
 RA_image -  RA_ref = -0.300"±0.013"
DEC_image - DEC_ref = -0.047"±0.008"
     rotation theta = 0.093±0.003 deg
 18 >> jw01309025001_02101_00010_nrca 
 RA_image -  RA_ref = -0.412"±0.012"
DEC_image - DEC_ref = 0.227"±0.010"
     rotation theta = 0.102±0.003 deg
 19 >> jw01309025001_02101_00010_nrcb 
 RA_image -  RA_ref = -0.296"±0.010"
DEC_image - DEC_ref = -0.047"±0.009"
     rotation theta = 0.092±0.003 deg
 20 >> jw01309025001_02101_00011_nrca 
 RA_image -  RA_ref = -0.417"±0.012"
DEC_image - DEC_ref = 0.227"±0.008"
     rotation theta = 0.098±0.003 deg
 21 >> jw01309025001_02101_00011_nrcb 
 RA_image -  RA_ref = -0.297"±0.013"
DEC_image - DEC_ref = -0.045"±0.009"
     rotation theta = 0.095±0.003 deg
 22 >> jw01309025001_02101_00012_nrca 
 RA_image -  RA_ref = -0.418"±0.012"
DEC_image - DEC_ref = 0.232"±0.008"
     rotation theta = 0.102±0.003 deg
 23 >> jw01309025001_02101_00012_nrcb 
 RA_image -  RA_ref = -0.297"±0.010"
DEC_image - DEC_ref = -0.043"±0.009"
     rotation theta = 0.095±0.003 deg
 24 >> jw01309025002_02101_00001_nrca 
 RA_image -  RA_ref = -0.411"±0.009"
DEC_image - DEC_ref = 0.295"±0.009"
     rotation theta = 0.097±0.002 deg
 25 >> jw01309025002_02101_00001_nrcb 
 RA_image -  RA_ref = -0.292"±0.014"
DEC_image - DEC_ref = 0.018"±0.008"
     rotation theta = 0.091±0.003 deg
 26 >> jw01309025002_02101_00002_nrca 
 RA_image -  RA_ref = -0.413"±0.008"
DEC_image - DEC_ref = 0.292"±0.009"
     rotation theta = 0.095±0.002 deg
 27 >> jw01309025002_02101_00002_nrcb 
 RA_image -  RA_ref = -0.293"±0.014"
DEC_image - DEC_ref = 0.015"±0.008"
     rotation theta = 0.090±0.003 deg
 28 >> jw01309025002_02101_00003_nrca 
 RA_image -  RA_ref = -0.412"±0.008"
DEC_image - DEC_ref = 0.290"±0.009"
     rotation theta = 0.096±0.002 deg
 29 >> jw01309025002_02101_00003_nrcb 
 RA_image -  RA_ref = -0.293"±0.014"
DEC_image - DEC_ref = 0.014"±0.008"
     rotation theta = 0.091±0.003 deg
 30 >> jw01309025002_02101_00004_nrca 
 RA_image -  RA_ref = -0.411"±0.008"
DEC_image - DEC_ref = 0.290"±0.008"
     rotation theta = 0.096±0.002 deg
 31 >> jw01309025002_02101_00004_nrcb 
 RA_image -  RA_ref = -0.293"±0.014"
DEC_image - DEC_ref = 0.013"±0.009"
     rotation theta = 0.092±0.003 deg
 32 >> jw01309025002_02101_00005_nrca 
 RA_image -  RA_ref = -0.430"±0.009"
DEC_image - DEC_ref = 0.296"±0.008"
     rotation theta = 0.094±0.002 deg
 33 >> jw01309025002_02101_00005_nrcb 
 RA_image -  RA_ref = -0.308"±0.012"
DEC_image - DEC_ref = 0.020"±0.008"
     rotation theta = 0.089±0.003 deg
 34 >> jw01309025002_02101_00006_nrca 
 RA_image -  RA_ref = -0.430"±0.008"
DEC_image - DEC_ref = 0.303"±0.008"
     rotation theta = 0.092±0.002 deg
 35 >> jw01309025002_02101_00006_nrcb 
 RA_image -  RA_ref = -0.306"±0.012"
DEC_image - DEC_ref = 0.022"±0.008"
     rotation theta = 0.090±0.003 deg
 36 >> jw01309025002_02101_00007_nrca 
 RA_image -  RA_ref = -0.427"±0.009"
DEC_image - DEC_ref = 0.303"±0.008"
     rotation theta = 0.094±0.002 deg
 37 >> jw01309025002_02101_00007_nrcb 
 RA_image -  RA_ref = -0.304"±0.012"
DEC_image - DEC_ref = 0.023"±0.008"
     rotation theta = 0.091±0.003 deg
 38 >> jw01309025002_02101_00008_nrca 
 RA_image -  RA_ref = -0.432"±0.008"
DEC_image - DEC_ref = 0.293"±0.008"
     rotation theta = 0.094±0.002 deg
 39 >> jw01309025002_02101_00008_nrcb 
 RA_image -  RA_ref = -0.308"±0.013"
DEC_image - DEC_ref = 0.017"±0.008"
     rotation theta = 0.085±0.002 deg
 40 >> jw01309025002_02101_00009_nrca 
 RA_image -  RA_ref = -0.443"±0.009"
DEC_image - DEC_ref = 0.300"±0.008"
     rotation theta = 0.097±0.002 deg
 41 >> jw01309025002_02101_00009_nrcb 
 RA_image -  RA_ref = -0.317"±0.011"
DEC_image - DEC_ref = 0.025"±0.008"
     rotation theta = 0.088±0.002 deg
 42 >> jw01309025002_02101_00010_nrca 
 RA_image -  RA_ref = -0.444"±0.010"
DEC_image - DEC_ref = 0.300"±0.008"
     rotation theta = 0.098±0.002 deg
 43 >> jw01309025002_02101_00010_nrcb 
 RA_image -  RA_ref = -0.318"±0.014"
DEC_image - DEC_ref = 0.029"±0.009"
     rotation theta = 0.087±0.002 deg
 44 >> jw01309025002_02101_00011_nrca 
 RA_image -  RA_ref = -0.444"±0.009"
DEC_image - DEC_ref = 0.302"±0.008"
     rotation theta = 0.096±0.002 deg
 45 >> jw01309025002_02101_00011_nrcb 
 RA_image -  RA_ref = -0.319"±0.012"
DEC_image - DEC_ref = 0.026"±0.009"
     rotation theta = 0.089±0.003 deg
 46 >> jw01309025002_02101_00012_nrca 
 RA_image -  RA_ref = -0.445"±0.009"
DEC_image - DEC_ref = 0.297"±0.009"
     rotation theta = 0.098±0.002 deg
 47 >> jw01309025002_02101_00012_nrcb 
 RA_image -  RA_ref = -0.321"±0.010"
DEC_image - DEC_ref = 0.022"±0.009"
     rotation theta = 0.087±0.002 deg
 48 >> jw01309026001_02101_00001_nrca 
 RA_image -  RA_ref = 0.265"±0.013"
DEC_image - DEC_ref = -0.067"±0.010"
     rotation theta = -0.115±0.003 deg
 49 >> jw01309026001_02101_00001_nrcb 
 RA_image -  RA_ref = 0.153"±0.012"
DEC_image - DEC_ref = 0.276"±0.009"
     rotation theta = -0.110±0.003 deg
 50 >> jw01309026001_02101_00002_nrca 
 RA_image -  RA_ref = 0.266"±0.013"
DEC_image - DEC_ref = -0.073"±0.009"
     rotation theta = -0.116±0.003 deg
 51 >> jw01309026001_02101_00002_nrcb 
 RA_image -  RA_ref = 0.155"±0.014"
DEC_image - DEC_ref = 0.270"±0.009"
     rotation theta = -0.114±0.004 deg
 52 >> jw01309026001_02101_00003_nrca 
 RA_image -  RA_ref = 0.265"±0.013"
DEC_image - DEC_ref = -0.064"±0.010"
     rotation theta = -0.115±0.003 deg
 53 >> jw01309026001_02101_00003_nrcb 
 RA_image -  RA_ref = 0.152"±0.012"
DEC_image - DEC_ref = 0.275"±0.009"
     rotation theta = -0.109±0.003 deg
 54 >> jw01309026001_02101_00004_nrca 
 RA_image -  RA_ref = 0.266"±0.014"
DEC_image - DEC_ref = -0.067"±0.010"
     rotation theta = -0.116±0.003 deg
 55 >> jw01309026001_02101_00004_nrcb 
 RA_image -  RA_ref = 0.155"±0.012"
DEC_image - DEC_ref = 0.275"±0.009"
     rotation theta = -0.111±0.003 deg
 56 >> jw01309026001_02101_00005_nrca 
 RA_image -  RA_ref = 0.287"±0.013"
DEC_image - DEC_ref = -0.071"±0.009"
     rotation theta = -0.112±0.002 deg
 57 >> jw01309026001_02101_00005_nrcb 
 RA_image -  RA_ref = 0.175"±0.016"
DEC_image - DEC_ref = 0.263"±0.012"
     rotation theta = -0.113±0.003 deg
 58 >> jw01309026001_02101_00006_nrca 
 RA_image -  RA_ref = 0.289"±0.015"
DEC_image - DEC_ref = -0.074"±0.009"
     rotation theta = -0.114±0.003 deg
 59 >> jw01309026001_02101_00006_nrcb 
 RA_image -  RA_ref = 0.183"±0.013"
DEC_image - DEC_ref = 0.266"±0.012"
     rotation theta = -0.115±0.003 deg
 60 >> jw01309026001_02101_00007_nrca 
 RA_image -  RA_ref = 0.288"±0.014"
DEC_image - DEC_ref = -0.075"±0.009"
     rotation theta = -0.115±0.003 deg
 61 >> jw01309026001_02101_00007_nrcb 
 RA_image -  RA_ref = 0.174"±0.013"
DEC_image - DEC_ref = 0.267"±0.011"
     rotation theta = -0.115±0.003 deg
 62 >> jw01309026001_02101_00008_nrca 
 RA_image -  RA_ref = 0.289"±0.013"
DEC_image - DEC_ref = -0.074"±0.009"
     rotation theta = -0.113±0.003 deg
 63 >> jw01309026001_02101_00008_nrcb 
 RA_image -  RA_ref = 0.174"±0.013"
DEC_image - DEC_ref = 0.267"±0.010"
     rotation theta = -0.114±0.003 deg
 64 >> jw01309026001_02101_00009_nrca 
 RA_image -  RA_ref = 0.314"±0.015"
DEC_image - DEC_ref = -0.083"±0.008"
     rotation theta = -0.109±0.003 deg
 65 >> jw01309026001_02101_00009_nrcb 
 RA_image -  RA_ref = 0.206"±0.011"
DEC_image - DEC_ref = 0.257"±0.009"
     rotation theta = -0.103±0.002 deg
 66 >> jw01309026001_02101_00010_nrca 
 RA_image -  RA_ref = 0.315"±0.013"
DEC_image - DEC_ref = -0.084"±0.009"
     rotation theta = -0.110±0.003 deg
 67 >> jw01309026001_02101_00010_nrcb 
 RA_image -  RA_ref = 0.205"±0.011"
DEC_image - DEC_ref = 0.254"±0.009"
     rotation theta = -0.106±0.002 deg
 68 >> jw01309026001_02101_00011_nrca 
 RA_image -  RA_ref = 0.310"±0.016"
DEC_image - DEC_ref = -0.083"±0.008"
     rotation theta = -0.111±0.003 deg
 69 >> jw01309026001_02101_00011_nrcb 
 RA_image -  RA_ref = 0.201"±0.012"
DEC_image - DEC_ref = 0.256"±0.012"
     rotation theta = -0.115±0.003 deg
 70 >> jw01309026001_02101_00012_nrca 
 RA_image -  RA_ref = 0.317"±0.013"
DEC_image - DEC_ref = -0.083"±0.008"
     rotation theta = -0.110±0.003 deg
 71 >> jw01309026001_02101_00012_nrcb 
 RA_image -  RA_ref = 0.201"±0.011"
DEC_image - DEC_ref = 0.257"±0.010"
     rotation theta = -0.110±0.003 deg
 72 >> jw01309026001_02106_00001_nrca 
 RA_image -  RA_ref = 0.265"±0.013"
DEC_image - DEC_ref = -0.066"±0.010"
     rotation theta = -0.115±0.003 deg
 73 >> jw01309026001_02106_00001_nrcb 
 RA_image -  RA_ref = 0.153"±0.012"
DEC_image - DEC_ref = 0.276"±0.009"
     rotation theta = -0.110±0.003 deg
 74 >> jw01309026001_02106_00002_nrca 
 RA_image -  RA_ref = 0.265"±0.013"
DEC_image - DEC_ref = -0.073"±0.009"
     rotation theta = -0.116±0.002 deg
 75 >> jw01309026001_02106_00002_nrcb 
 RA_image -  RA_ref = 0.155"±0.014"
DEC_image - DEC_ref = 0.272"±0.009"
     rotation theta = -0.113±0.003 deg
 76 >> jw01309026001_02106_00003_nrca 
 RA_image -  RA_ref = 0.265"±0.013"
DEC_image - DEC_ref = -0.064"±0.010"
     rotation theta = -0.115±0.003 deg
 77 >> jw01309026001_02106_00003_nrcb 
 RA_image -  RA_ref = 0.152"±0.013"
DEC_image - DEC_ref = 0.276"±0.009"
     rotation theta = -0.109±0.003 deg
 78 >> jw01309026001_02106_00004_nrca 
 RA_image -  RA_ref = 0.266"±0.013"
DEC_image - DEC_ref = -0.068"±0.010"
     rotation theta = -0.115±0.003 deg
 79 >> jw01309026001_02106_00004_nrcb 
 RA_image -  RA_ref = 0.155"±0.012"
DEC_image - DEC_ref = 0.275"±0.009"
     rotation theta = -0.110±0.003 deg
 80 >> jw01309026001_02106_00005_nrca 
 RA_image -  RA_ref = 0.286"±0.014"
DEC_image - DEC_ref = -0.075"±0.009"
     rotation theta = -0.113±0.002 deg
 81 >> jw01309026001_02106_00005_nrcb 
 RA_image -  RA_ref = 0.175"±0.014"
DEC_image - DEC_ref = 0.263"±0.012"
     rotation theta = -0.108±0.003 deg
 82 >> jw01309026001_02106_00006_nrca 
 RA_image -  RA_ref = 0.289"±0.015"
DEC_image - DEC_ref = -0.073"±0.009"
     rotation theta = -0.113±0.003 deg
 83 >> jw01309026001_02106_00006_nrcb 
 RA_image -  RA_ref = 0.178"±0.013"
DEC_image - DEC_ref = 0.265"±0.012"
     rotation theta = -0.115±0.003 deg
 84 >> jw01309026001_02106_00007_nrca 
 RA_image -  RA_ref = 0.288"±0.015"
DEC_image - DEC_ref = -0.071"±0.009"
     rotation theta = -0.113±0.002 deg
 85 >> jw01309026001_02106_00007_nrcb 
 RA_image -  RA_ref = 0.174"±0.013"
DEC_image - DEC_ref = 0.266"±0.010"
     rotation theta = -0.115±0.003 deg
 86 >> jw01309026001_02106_00008_nrca 
 RA_image -  RA_ref = 0.289"±0.012"
DEC_image - DEC_ref = -0.074"±0.009"
     rotation theta = -0.114±0.003 deg
 87 >> jw01309026001_02106_00008_nrcb 
 RA_image -  RA_ref = 0.174"±0.013"
DEC_image - DEC_ref = 0.266"±0.010"
     rotation theta = -0.106±0.003 deg
 88 >> jw01309026001_02106_00009_nrca 
 RA_image -  RA_ref = 0.314"±0.015"
DEC_image - DEC_ref = -0.083"±0.008"
     rotation theta = -0.110±0.003 deg
 89 >> jw01309026001_02106_00009_nrcb 
 RA_image -  RA_ref = 0.207"±0.011"
DEC_image - DEC_ref = 0.257"±0.009"
     rotation theta = -0.104±0.002 deg
 90 >> jw01309026001_02106_00010_nrca 
 RA_image -  RA_ref = 0.316"±0.013"
DEC_image - DEC_ref = -0.086"±0.009"
     rotation theta = -0.112±0.003 deg
 91 >> jw01309026001_02106_00010_nrcb 
 RA_image -  RA_ref = 0.206"±0.011"
DEC_image - DEC_ref = 0.255"±0.010"
     rotation theta = -0.107±0.002 deg
 92 >> jw01309026001_02106_00011_nrca 
 RA_image -  RA_ref = 0.310"±0.015"
DEC_image - DEC_ref = -0.083"±0.009"
     rotation theta = -0.111±0.003 deg
 93 >> jw01309026001_02106_00011_nrcb 
 RA_image -  RA_ref = 0.198"±0.011"
DEC_image - DEC_ref = 0.256"±0.010"
     rotation theta = -0.115±0.003 deg
 94 >> jw01309026001_02106_00012_nrca 
 RA_image -  RA_ref = 0.318"±0.012"
DEC_image - DEC_ref = -0.083"±0.008"
     rotation theta = -0.110±0.003 deg
 95 >> jw01309026001_02106_00012_nrcb 
 RA_image -  RA_ref = 0.201"±0.011"
DEC_image - DEC_ref = 0.257"±0.010"
     rotation theta = -0.110±0.003 deg
In [ ]:
 
In [40]:
'''Checking the derived rotation, if necessary'''
fig, ax = plt.subplots(1, 2, figsize = (12, 6) )
ax[0].plot(tmp_coord_daofind[idx_daofind].dec, tmp_ra_offset, marker = '.', color = 'k', ls = 'none')
ax[0].set_ylim(-0.4, 0.4)


ax[1].plot(tmp_coord_daofind[idx_daofind].ra, tmp_dec_offset, marker = '.', color = 'k', ls = 'none')
ax[1].set_ylim(-0.4, 0.4)


poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec, ydata = tmp_ra_offset)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra, ydata = tmp_dec_offset)[0]
## iter-1
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])[0]
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])[0]
## iter-2
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA, pcov_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC, pccov_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])
perr_dec_dRA = np.sqrt(np.diag(pcov_dec_dRA))
perr_ra_dDEC = np.sqrt(np.diag(pccov_ra_dDEC))

ax[0].plot(tmp_coord_daofind[idx_daofind].dec.value, 
           np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value))
ax[1].plot(tmp_coord_daofind[idx_daofind].ra.value, 
           np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value))
ax[0].set(xlabel = 'Dec', ylabel = 'dRA')
ax[1].set(xlabel = 'RA', ylabel = 'dDEC')
plt.tight_layout()

print(poly_dec_dRA, poly_ra_dDEC)
[  -7.26053243 -561.80868745] [   1.42702112 -237.44536921]

Write astrometric catalog if necessary. The one that I got for Ice Age data (F200W image & F444W grism) is available at http://gxn.as.arizona.edu/~sunfengwu/IceAge/PID_1309_SW_wavecal_F200W_F430M.dat.

In [43]:
# ascii.write(tb_sw_astrometry, '/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F200W_F430M.dat', 
#             overwrite = True)
In [42]:
# tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_relative_astrom_obs24.dat')

print('0-16-50-84-100 percentiles of astrometric accuracy for RA direction [in mas]:')
print(np.round(np.percentile(tb_sw_astrometry['dRA_err'].data / np.sqrt(tb_sw_astrometry['N_match'].data - 1), 
                             [0, 16, 50, 84, 100]) * 1e3, 2))
0-16-50-84-100 percentiles of astrometric accuracy for RA direction [in mas]:
[ 8.56 10.34 13.93 15.99 19.05]
In [46]:
'''Evolution of astrometric offsets with exposures:'''
fig = plt.figure(constrained_layout = False, figsize = (12, 6))
ax_dict = fig.subplot_mosaic('''AACC
                                BBCC''')
ax_dict['A'].plot(tb_sw_astrometry['dRA']) # color =  plt.cm.Reds(0.5))
ax_dict['B'].plot(tb_sw_astrometry['dDEC']) # color =  plt.cm.Blues(0.5))
ax_dict['A'].set(ylabel = 'dRA [arcsec]', ylim = (-0.6, 0.6))
ax_dict['B'].set(ylabel = 'dDec [arcsec]', ylim = (-0.6, 0.6), xlabel = '$i_\mathrm{frame}$')

for i, tmp_dRA in enumerate(tb_sw_astrometry['dRA']):
    if i%2 == 0: # mod-A
        kwargs_dot = dict(color = plt.cm.Reds(i/len(tb_sw_astrometry)), marker = 'o', zorder = i)
    else:        # mod-B
        kwargs_dot = dict(color = plt.cm.Blues(i/len(tb_sw_astrometry)), marker = 's', zorder = i)
    ax_dict['C'].plot(tb_sw_astrometry['dRA'][i], tb_sw_astrometry['dDEC'][i], **kwargs_dot, 
                      mec = 'silver', mew = 1, ms = 8)
        
ax_dict['C'].plot(tb_sw_astrometry['dRA'][::2], tb_sw_astrometry['dDEC'][::2], 
                  ls = '-', color = plt.cm.Reds(0.5), zorder = 1)
ax_dict['C'].plot(tb_sw_astrometry['dRA'][1::2], tb_sw_astrometry['dDEC'][1::2], 
                  ls = '-', color = plt.cm.Blues(0.5), zorder = 0)
ax_dict['C'].set(aspect = 1, ylabel = 'dDEC [arcsec]', xlabel = 'dRA [arcsec]')
corner_text(ax_dict['C'], s = 'modA', color =  plt.cm.Reds(0.5), loc = 2, weight = 'semibold')
corner_text(ax_dict['C'], s = 'modB', color =  plt.cm.Blues(0.5), loc = 1, weight = 'semibold')
fig.text(0.5, 0.98, 'PID-1309 SW Relative Astrometric Offset', ha = 'center', va = 'top', fontsize = 18)
plt.tight_layout()
plt.subplots_adjust(top = 0.93)
# fig.savefig('/home/u24/fengwusun/IceAge/PID_1309_SW_relative_astrom_F200W.pdf', dpi = 100)
In [47]:
'''Astrometric offsets between module A and B'''
plt.errorbar(np.arange(len(tb_sw_astrometry)/2),
             ((tb_sw_astrometry['dRA'][::2] - tb_sw_astrometry['dRA'][1::2])**2 + 
              (tb_sw_astrometry['dDEC'][::2] - tb_sw_astrometry['dDEC'][1::2])**2)**0.5,
             yerr = tb_sw_astrometry['dDEC_err'][::2] * np.sqrt(2) / np.sqrt(tb_sw_astrometry['N_match'][::2] - 1),
         ls = '-', zorder = 0)
plt.gca().set(xlabel = 'Exp NUM', ylabel = 'modA - modB OFFSET', ylim = (-0.2, 0.5))
Out[47]:
[Text(0.5, 0, 'Exp NUM'), Text(0, 0.5, 'modA - modB OFFSET'), (-0.2, 0.5)]

Check source position predicted by fits WCS + Measured Astrometry on SW images

This method will also be used to predict source positions on grism data

In [50]:
## Select one SW exposure to look at:
k = 0
tmp_img_cal_path = list_img_cal_this_band[k * 4 + 1]
tmp_img_cal_path_base = os.path.basename(tmp_img_cal_path)[:30]
## astrometric info:
item_sw_astrom = tb_sw_astrometry[tb_sw_astrometry['expName'] == tmp_img_cal_path_base][0]
print('%3d >> %s ' % (k, tmp_img_cal_path))
tmp_img_cal_hd = fits.getheader(tmp_img_cal_path, 0)
tmp_detector = tmp_img_cal_hd['DETECTOR'].lower()
tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
tmp_img_cal_wcs = wcs.WCS(tmp_img_cal_sci_hd)
time_obs = Time(tmp_img_cal_sci_hd['MJD-AVG'], format ='mjd')
tmp_detectors = np.array(['%s%d' % (tmp_detector[:-1], tmp_det) for tmp_det in [1, 2, 3, 4]])
tmp_img_cal_paths = np.array([tmp_img_cal_path.replace(tmp_detector, x) for x in tmp_detectors])
    
'''astrometric catalog loaded'''
# tmp_RA, tmp_DEC = tb_LW['ALPHA_J2000'].data, tb_LW['DELTA_J2000'].data 
tmp_RA, tmp_DEC = tb_LW['RA'].data, tb_LW['DEC'].data 
# tmp_RA_astrom = tmp_RA + item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_DEC)) / 3600.
# tmp_DEC_astrom = tmp_DEC + item_sw_astrom['dDEC'] / 3600.
# tmp_img_cal_sci_hd['crval1'] -= item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_img_cal_sci_hd['crval2'])) / 3600.
# tmp_img_cal_sci_hd['crval2'] -= item_sw_astrom['dDEC'] / 3600.
'''astrometric info applied on WCS:'''
tmp_cd_matrix = np.array([[tmp_img_cal_sci_hd['CD1_1'], tmp_img_cal_sci_hd['CD1_2']],
                          [tmp_img_cal_sci_hd['CD2_1'], tmp_img_cal_sci_hd['CD2_2']]])
tmp_rot_matrix = np.array([[np.cos(np.deg2rad(-item_sw_astrom['theta'])), - np.sin(np.deg2rad(-item_sw_astrom['theta']))],
                           [np.sin(np.deg2rad(-item_sw_astrom['theta'])),   np.cos(np.deg2rad(-item_sw_astrom['theta']))]])
tmp_cd_matrix = np.matmul(tmp_cd_matrix, tmp_rot_matrix)
tmp_img_cal_sci_hd['CD1_1'] = tmp_cd_matrix[0][0]
tmp_img_cal_sci_hd['CD1_2'] = tmp_cd_matrix[0][1]
tmp_img_cal_sci_hd['CD2_1'] = tmp_cd_matrix[1][0]
tmp_img_cal_sci_hd['CD2_2'] = tmp_cd_matrix[1][1]

                          
# ## Read latest SiAF for NIRCam
# nircam_siaf = pysiaf.Siaf('NIRCam')
# tmp_siaf = nircam_siaf['%s_%s' % (tmp_img_cal_hd['detector'], tmp_img_cal_hd['SUBARRAY'])]
# try: rotation = item_sw_astrom['theta']
# except KeyError: rotation = 0
    
'''
Convert coordinate of astrometric sources to pixel_x/y based on WCS
'''
plt.figure(figsize = (12, 12))
tmp_sci_img = fits.getdata(tmp_img_cal_path, 'sci')
plt.imshow(tmp_sci_img, vmin = np.nanpercentile(tmp_sci_img, 40), vmax = np.nanpercentile(tmp_sci_img, 98), origin = 'lower')


tmp_wcs = wcs.WCS(tmp_img_cal_sci_hd)
pixelx, pixely = wcs.utils.skycoord_to_pixel(SkyCoord(tmp_RA, tmp_DEC, unit = (u.deg, u.deg)), tmp_wcs)  
    
arg_gaia_plot = np.where((pixelx > 0) & (pixelx < 2048) & (pixely > 0) & (pixely < 2048))[0]
plt.plot(pixelx[arg_gaia_plot], pixely[arg_gaia_plot], marker = 'o', mfc = 'none', mec = 'salmon', ms = 10, ls = 'none', mew = 2)


# tmp_asdf_distortion = asdf.open(crds.getreferences(tmp_img_cal_hd, 
#                                                    reftypes = ['distortion'], ignore_cache = False)['distortion'])
# loc_v2, loc_v3 = tmp_asdf_distortion['model'](pixelx, pixely)
# pixelx, pixely = tmp_siaf.tel_to_sci(loc_v2, loc_v3) 
# plt.plot(pixelx[arg_gaia_plot], pixely[arg_gaia_plot], marker = 'o', mfc = 'none', mec = 'w', ms = 10, ls = 'none', mew = 2)

# plt.xlim(200, 700);plt.ylim(200, 700)
# plt.xlim(1300, 1800);plt.ylim(1300, 1800)
  0 >> /xdisk/egami/fengwusun/ERSs/IceAge/grism/../direct/jw01309025001_02101_00001_nrca2_cal.fits 
Out[50]:
[<matplotlib.lines.Line2D at 0x7f83354d6310>]
In [ ]:
 

5. Genertae Source Catalog for each LW grism exposure

With astrometric information in hand, we can start to generate the catalog that include sources of our interest for each LW grism exposures. With these catalogs (I name them as "POM-applied catalogs" because POM selection effects are applied), you can then run the code blocks to extract 2D spectra from individual LW grism frames, co-add them together to produce the final 2D grism spectra.

First of all, let us grab the reduced LW grism images (_rate_lv1.5.fits) again:

In [51]:
'''LW Grism images and corresponding astrometric information:'''
# all_rate_list = glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[3-4]*_rate_lv1.5.fits')
# tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F150W_F410M.dat')
all_rate_list = glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[5-6]*_rate_lv1.5.fits')
tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F200W_F430M.dat')

all_rate_list.sort()
all_rate_list = np.array(all_rate_list)
print(len(all_rate_list), 'wcs corrected fits files')

'''read filter, module and pupil'''
all_filters, all_module, all_pupil = [], [], []
all_roll_angle, all_target = [], []
for i, x in enumerate(all_rate_list):
    tmp_header = fits.getheader(x)
    all_filters.append(tmp_header['filter'])
    all_module.append(tmp_header['module'])
    all_pupil.append(tmp_header['pupil'])
    all_roll_angle.append(fits.getheader(x, 1)['ROLL_REF'])
    all_target.append(tmp_header['TARGPROP'])
all_filters, all_module, all_pupil = np.array(all_filters), np.array(all_module), np.array(all_pupil)
all_roll_angle, all_target = np.array(all_roll_angle), np.array(all_target)
print(np.unique(all_filters), np.unique(all_module), np.unique(all_target))

### change target prefix if you like:
tmp_filter = np.unique(all_filters)[0]
if tmp_filter == 'F322W2': tmp_target = 'CHAMMS1-H2O'
if tmp_filter == 'F444W': tmp_target = 'CHAMMS1-CO'
96 wcs corrected fits files
['F444W'] ['A' 'B'] ['CHAMMS1-C2-FIELD']

Ice Age data consist of F322W2 and F444W grism spectra taken with both dispersion direction (R and C). In the following example, we input the source catalog extracted in both F410M (imaged after F322W2 grism exposure) and F430M (imaged after F444W grism exposure) band as master catalog for spectral extraction.

In [53]:
'''Input master catalog for spectral extraction:'''
tb_f322w2 = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F410M_source.cat')
tb_f444w = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F430M_source.cat')
tb_f322w2 = tb_f322w2[tb_f322w2['flags'] == 0]
tb_f444w = tb_f444w[tb_f444w['flags'] == 0]
'''If necessary, select source based on morphology and brigtness'''
# tb_f322w2 = tb_f322w2[(tb_f322w2['MAG_AUTO'] < 18.5) & (tb_f322w2['CLASS_STAR'] > 0.5)]
### read coordniates
coord_f322w2 = SkyCoord(tb_f322w2['RA'], tb_f322w2['DEC'], unit = (u.deg, u.deg))
coord_f444w = SkyCoord(tb_f444w['RA'], tb_f444w['DEC'], unit = (u.deg, u.deg))
print(len(coord_f322w2), len(coord_f444w))
1032 649
In [55]:
'''Columnames of my catalog'''
print(tb_f444w.colnames)
tb_f444w[:4]['id'].data
['id', 'xcentroid', 'ycentroid', 'sharpness', 'roundness1', 'roundness2', 'npix', 'sky', 'peak', 'flux', 'mag', 'RA', 'DEC', 'flux_r2', 'flux_r4', 'fluxerr_r4', 'dao_indicator', 'aper_indicator', 'SN_r4', 'flags', 'flux_r4_corr', 'fluxerr_r4_corr', 'abmag', 'abmag_err']
Out[55]:
array([ 521,  498,  805, 1540])

For each exposure, apply astrometric correction and generate source catalog for spectral extraction

In this step, there is one tricky thing called filter offset. As you might have recalled, these offsets are introduced by different SW+LW filter combination. In commissioning, we used F212N + F444W for F444W grism spectroscopy, but Ice Age used F200W + F444W. Therefore, the spectral tracing model derived in commissioning could be offset from the actual ones in IceAge data due to the difference in SW filters. This should be a global shift (dx, dy) which is different for different modules but consistent for all exposures. I believe that this has not been well calibrated so far, so you will need to enter some dx, dy value manually. You can use the script a few blocks below the following code block to check the filter offsets in current CRDS, but my general advice is to change these dx, dy number manually after you successfully extract some 2D spectra and measure the real offsets.

In [60]:
'''POM applied catalogs'''
### `all_rate_list`: a list of paths to the grism exposure fits files
for i, tmp_rate_path in enumerate(all_rate_list[:]):
    ''' 
    ## The following commands are used to read the astrometric offsets which is measured with SW detectors.
    ## In real life, the pointing accuracy of the telescope may not be stable and the absolute RA and DEC 
    ## offset can be floating. We will need to measure the astrometric shift from the SW exposures associated 
    ## with the grism exposures, and correct this factor for the RA and DEC in each source catalog associated 
    ## with the grism exposures. 
    '''
    tmp_rate_path_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
    print('[%3d]' % i, tmp_rate_path_base)
    
    tmp_grism_hd_1st = fits.getheader(tmp_rate_path, 0)
    tmp_grism_hd_sci = fits.getheader(tmp_rate_path, 'sci')
    
    tmp_filter = tmp_grism_hd_1st['filter']
    tmp_loc_v2, tmp_loc_v3 = tmp_grism_hd_sci['V2_REF'], tmp_grism_hd_sci['V3_REF']
    ### for different grism filters, read different catalog
    if tmp_filter == 'F444W': list_coords, tb_sex = coord_f444w, tb_f444w
    elif tmp_filter == 'F322W2': list_coords, tb_sex = coord_f322w2, tb_f322w2
        
    # if tmp_grism_hd_1st['filter'] != 'F444W': continue
    
    '''Accurate RA & Dec of Direct-Imaging sources  '''
    tmp_RA, tmp_DEC = list_coords.ra.to(u.deg).value, list_coords.dec.to(u.deg).value
    ### if the catalog is gaia star, considering proper motion:
    # time_obs = Time(tmp_grism_hd_sci['MJD-AVG'], format ='mjd')
    # dRA = ((time_obs - time_gaia).to(u.yr).value * tb_gaia['pmRA'].data * u.mas / np.cos(np.deg2rad(tb_gaia['DE_ICRS']))).to(u.deg).value
    # dDec = ((time_obs - time_gaia).to(u.yr).value * tb_gaia['pmDE'].data * u.mas).to(u.deg).value
    # tmp_RA, tmp_DEC = tb_gaia['RA_ICRS'].data.data + dRA, tb_gaia['DE_ICRS'].data.data + dDec
    '''Apply astrometric error measured on SW images'''
    item_sw_astrom = tb_sw_astrometry[tb_sw_astrometry['expName'] == tmp_rate_path_base]
    if len(item_sw_astrom) > 0:
        ## if there is astrometric information, consider that
        item_sw_astrom = item_sw_astrom[0]
        tmp_grism_hd_sci['CRVAL1'] -= item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_grism_hd_sci['CRVAL2'])) / 3600.
        tmp_grism_hd_sci['CRVAL2'] -= item_sw_astrom['dDEC'] / 3600. 
        tmp_cd_matrix = np.array([[tmp_grism_hd_sci['CD1_1'], tmp_grism_hd_sci['CD1_2']],
                                  [tmp_grism_hd_sci['CD2_1'], tmp_grism_hd_sci['CD2_2']]])
        tmp_rot_matrix = np.array([[np.cos(np.deg2rad(-item_sw_astrom['theta'])), - np.sin(np.deg2rad(-item_sw_astrom['theta']))],
                                   [np.sin(np.deg2rad(-item_sw_astrom['theta'])),   np.cos(np.deg2rad(-item_sw_astrom['theta']))]])
        tmp_cd_matrix = np.matmul(tmp_cd_matrix, tmp_rot_matrix)
        tmp_grism_hd_sci['CD1_1'] = tmp_cd_matrix[0][0]
        tmp_grism_hd_sci['CD1_2'] = tmp_cd_matrix[0][1]
        tmp_grism_hd_sci['CD2_1'] = tmp_cd_matrix[1][0]
        tmp_grism_hd_sci['CD2_2'] = tmp_cd_matrix[1][1]
        tmp_RA_astrom = tmp_RA # + item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_DEC)) / 3600.
        tmp_DEC_astrom = tmp_DEC # + item_sw_astrom['dDEC'] / 3600.
    else:
        ## if no astrometric information avaialble, skip this step and directly use RA and DEC
        print('Warning: No SW images found for %s, skip astrom correction! ' % tmp_rate_path_base)
        tmp_RA_astrom = tmp_RA 
        tmp_DEC_astrom = tmp_DEC 
    try: rotation = item_sw_astrom['theta']
    except KeyError: rotation = 0
    
    '''RA/DEC to grism-frame x and y'''
    tmp_grism_wcs = wcs.WCS(tmp_grism_hd_sci)
    idx_this_field = np.where((np.abs((tmp_RA_astrom - tmp_grism_hd_sci['crval1']) * np.cos(np.deg2rad(tmp_grism_hd_sci['crval2']))) < 3 / 60.) &
                              (np.abs(tmp_DEC_astrom - tmp_grism_hd_sci['crval2']) < 3 / 60.))[0]
    pixelx, pixely  = wcs.utils.skycoord_to_pixel(SkyCoord(tmp_RA_astrom[idx_this_field], 
                                                           tmp_DEC_astrom[idx_this_field], unit = (u.deg, u.deg)), 
                                                  tmp_grism_wcs)
    '''This is the place to add small offsets from all kind of reasons (e.g., filter offset)'''
    if tmp_grism_hd_1st['module'] == 'A' : dx, dy = 0, 0  ###  +2.0, -2.5
    else: dx, dy = 0, 0
    pixelx, pixely = pixelx + dx, pixely + dy
    tb_sex = tb_sex[idx_this_field]

    '''for sextractor-based catalog '''
#     tmp_mag_auto, tmp_magerr_auto = tb_sex['MAG_AUTO'].data, tb_sex['MAGERR_AUTO'].data
#     tb_pom_applied = Table(data = [tb_sex['NUMBER'].data, 
#                                    tmp_RA[idx_this_field], tmp_DEC[idx_this_field], 
#                                    pixelx, pixely, tmp_mag_auto, tmp_magerr_auto,
#                                    tb_sex['A_IMAGE'], tb_sex['B_IMAGE'],  tb_sex['THETA_IMAGE'], 
#                                    tb_sex['FLAGS'], tb_sex['CLASS_STAR']],
#                            names = ['Index', 'ra', 'dec', 'pixel_x', 'pixel_y', 'MAG_AUTO', 'MAGERR_AUTO',
#                                     'A_IMAGE', 'B_IMAGE', 'THETA_IMAGE', 'FLAGS', 'CLASS_STAR'])
    '''for DAOFIND-based catalog '''
    tmp_mag_auto, tmp_magerr_auto = tb_sex['abmag'].data, tb_sex['abmag_err'].data
    tb_pom_applied = Table(data = [tb_sex['id'].data, 
                                   tmp_RA[idx_this_field], tmp_DEC[idx_this_field], 
                                   pixelx, pixely, tmp_mag_auto, tmp_magerr_auto, tb_sex['flags']],
                           names = ['Index', 'ra', 'dec', 'pixel_x', 'pixel_y', 'MAG_AUTO', 'MAGERR_AUTO', 'FLAGS'])
    # arg_is_pickoff = np.where(is_pickoff(pixelx, pixely, module = tmp_grism_hd_1st['module']))[0]
    arg_is_pickoff = np.where(is_pickoff_PS(pixelx, pixely, module = tmp_grism_hd_1st['module'], 
                                            filter = tmp_grism_hd_1st['filter'], pupil = tmp_grism_hd_1st['pupil'][-1]) > 0.05)[0]
    tb_pom_applied = tb_pom_applied[arg_is_pickoff]
    tb_pom_applied['ra'].info.format = '.6f'
    tb_pom_applied['dec'].info.format = '.6f'
    tb_pom_applied['pixel_x'].info.format = '.3f'
    tb_pom_applied['pixel_y'].info.format = '.3f'
    tb_pom_applied['MAG_AUTO'].info.format = '.3f'
    tb_pom_applied['MAGERR_AUTO'].info.format = '.3f'
    try:
        tb_pom_applied['A_IMAGE'].info.format = '.3f'
        tb_pom_applied['B_IMAGE'].info.format = '.3f'
        tb_pom_applied['THETA_IMAGE'].info.format = '.2f'
        tb_pom_applied['CLASS_STAR'].info.format = '.3f'
    except KeyError: pass
    
    '''Change the path of catalogs to be saved:'''
    path_tb_pom_applied = '/home/u24/fengwusun/IceAge/POM_catalog/' +  tmp_rate_path_base + '_%s_dirimg_sources.list' % tmp_filter
    ascii.write(tb_pom_applied, path_tb_pom_applied, overwrite = True)
[  0] jw01309025001_02101_00001_nrca
[  1] jw01309025001_02101_00001_nrcb
[  2] jw01309025001_02101_00002_nrca
[  3] jw01309025001_02101_00002_nrcb
[  4] jw01309025001_02101_00003_nrca
[  5] jw01309025001_02101_00003_nrcb
[  6] jw01309025001_02101_00004_nrca
[  7] jw01309025001_02101_00004_nrcb
[  8] jw01309025001_02101_00005_nrca
[  9] jw01309025001_02101_00005_nrcb
[ 10] jw01309025001_02101_00006_nrca
[ 11] jw01309025001_02101_00006_nrcb
[ 12] jw01309025001_02101_00007_nrca
[ 13] jw01309025001_02101_00007_nrcb
[ 14] jw01309025001_02101_00008_nrca
[ 15] jw01309025001_02101_00008_nrcb
[ 16] jw01309025001_02101_00009_nrca
[ 17] jw01309025001_02101_00009_nrcb
[ 18] jw01309025001_02101_00010_nrca
[ 19] jw01309025001_02101_00010_nrcb
[ 20] jw01309025001_02101_00011_nrca
[ 21] jw01309025001_02101_00011_nrcb
[ 22] jw01309025001_02101_00012_nrca
[ 23] jw01309025001_02101_00012_nrcb
[ 24] jw01309025002_02101_00001_nrca
[ 25] jw01309025002_02101_00001_nrcb
[ 26] jw01309025002_02101_00002_nrca
[ 27] jw01309025002_02101_00002_nrcb
[ 28] jw01309025002_02101_00003_nrca
[ 29] jw01309025002_02101_00003_nrcb
[ 30] jw01309025002_02101_00004_nrca
[ 31] jw01309025002_02101_00004_nrcb
[ 32] jw01309025002_02101_00005_nrca
[ 33] jw01309025002_02101_00005_nrcb
[ 34] jw01309025002_02101_00006_nrca
[ 35] jw01309025002_02101_00006_nrcb
[ 36] jw01309025002_02101_00007_nrca
[ 37] jw01309025002_02101_00007_nrcb
[ 38] jw01309025002_02101_00008_nrca
[ 39] jw01309025002_02101_00008_nrcb
[ 40] jw01309025002_02101_00009_nrca
[ 41] jw01309025002_02101_00009_nrcb
[ 42] jw01309025002_02101_00010_nrca
[ 43] jw01309025002_02101_00010_nrcb
[ 44] jw01309025002_02101_00011_nrca
[ 45] jw01309025002_02101_00011_nrcb
[ 46] jw01309025002_02101_00012_nrca
[ 47] jw01309025002_02101_00012_nrcb
[ 48] jw01309026001_02101_00001_nrca
[ 49] jw01309026001_02101_00001_nrcb
[ 50] jw01309026001_02101_00002_nrca
[ 51] jw01309026001_02101_00002_nrcb
[ 52] jw01309026001_02101_00003_nrca
[ 53] jw01309026001_02101_00003_nrcb
[ 54] jw01309026001_02101_00004_nrca
[ 55] jw01309026001_02101_00004_nrcb
[ 56] jw01309026001_02101_00005_nrca
[ 57] jw01309026001_02101_00005_nrcb
[ 58] jw01309026001_02101_00006_nrca
[ 59] jw01309026001_02101_00006_nrcb
[ 60] jw01309026001_02101_00007_nrca
[ 61] jw01309026001_02101_00007_nrcb
[ 62] jw01309026001_02101_00008_nrca
[ 63] jw01309026001_02101_00008_nrcb
[ 64] jw01309026001_02101_00009_nrca
[ 65] jw01309026001_02101_00009_nrcb
[ 66] jw01309026001_02101_00010_nrca
[ 67] jw01309026001_02101_00010_nrcb
[ 68] jw01309026001_02101_00011_nrca
[ 69] jw01309026001_02101_00011_nrcb
[ 70] jw01309026001_02101_00012_nrca
[ 71] jw01309026001_02101_00012_nrcb
[ 72] jw01309026001_02106_00001_nrca
[ 73] jw01309026001_02106_00001_nrcb
[ 74] jw01309026001_02106_00002_nrca
[ 75] jw01309026001_02106_00002_nrcb
[ 76] jw01309026001_02106_00003_nrca
[ 77] jw01309026001_02106_00003_nrcb
[ 78] jw01309026001_02106_00004_nrca
[ 79] jw01309026001_02106_00004_nrcb
[ 80] jw01309026001_02106_00005_nrca
[ 81] jw01309026001_02106_00005_nrcb
[ 82] jw01309026001_02106_00006_nrca
[ 83] jw01309026001_02106_00006_nrcb
[ 84] jw01309026001_02106_00007_nrca
[ 85] jw01309026001_02106_00007_nrcb
[ 86] jw01309026001_02106_00008_nrca
[ 87] jw01309026001_02106_00008_nrcb
[ 88] jw01309026001_02106_00009_nrca
[ 89] jw01309026001_02106_00009_nrcb
[ 90] jw01309026001_02106_00010_nrca
[ 91] jw01309026001_02106_00010_nrcb
[ 92] jw01309026001_02106_00011_nrca
[ 93] jw01309026001_02106_00011_nrcb
[ 94] jw01309026001_02106_00012_nrca
[ 95] jw01309026001_02106_00012_nrcb

Optional: Check filter offset:

This should have been corrected in stage-2 using STScI CRDS files. However, I stil see some offsets in modA, weird!

In [57]:
### select a reduced grism exposure:
tmp_rate_path = '/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00012_nrcblong_rate_lv1.5.fits'
tmp_rate_path_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
### SW header
tmp_sw_cal_path = os.path.dirname(tmp_rate_path) + '/../direct/' + tmp_rate_path_base + '1_cal.fits'
tmp_sw_cal_hd = fits.getheader(tmp_sw_cal_path)
### LW header
tmp_grism_hd_1st = fits.getheader(tmp_rate_path, 0)
### offsets in LW & SW
tmp_tb_foff_lw = Table(asdf.open(crds.getreferences(tmp_grism_hd_1st)['filteroffset'])['filters'])
tmp_tb_foff_sw = Table(asdf.open(crds.getreferences(tmp_sw_cal_hd)['filteroffset'])['filters'])
### Filter offset of  F212N, which is used in Commissioning
item_F212N = tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == 'F212N') & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
### Filter offset of the filter used in the current SW exposures
item_sw = tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == tmp_sw_cal_hd['filter']) & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
In [59]:
item_F212N, item_sw
# tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == 'F200W') & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
Out[59]:
(<Row index=0>
 column_offset filter pupil row_offset
    float64     str6   str8  float64  
 ------------- ------ ----- ----------
          -0.1  F212N CLEAR      0.165,
 <Row index=0>
 column_offset filter pupil row_offset
    float64     str6   str8  float64  
 ------------- ------ ----- ----------
         1.016  F200W CLEAR      1.625)
In [739]:
'''Sometimes LW filters also introduce offset (e.g., F356W - F322W2). Check that if necessary'''
# tmp_tb_foff_lw[(tmp_tb_foff_lw['filter'] == 'F356W') & (tmp_tb_foff_lw['pupil'] == 'CLEAR')]

6. 2D Spectral Extraction

The following step is to extract grism spectra with the POM-applied catalogs and reduced grism data you produced above. At this step, you will need to load the following functions:

  1. Spectral Tracing Model

Now we assume $(x_0, y_0)$ are direct imaging source position, $(x_s, y_s)$ are spectral pixel position at wavelength $\lambda_s$, and $dx = x_s - x_0$, $dy = y_s - y_0$. The spectral tracing model gives $dy(x_0, y_0, dx)$, i.e., the spectral tracing. The following plot shows the spectral tracing that we derived for F322W2 module A grism R using the commissioning data taken for the LMC field. Red squares are direct imaging positions of sources taken in multiple exposures, and black points are the corresponding spectral traces (curvature is exaggerated). The plot of dx versus dy, color-coded by the source position, is shown on the right. Note the spectra with $dx \sim 1000$ are second-order, which have not been fully characterized at this step.

spectral tracing in F322W2 modA grismR

The spectral tracing function is defined as fit_disp_order23() (see the beginning of this notebook). The best-fit parameters for these functions in different filter/module/grism can be found at http://u.arizona.edu/~fengwusun/data/dy_dx_tracing_model/. For filters at 2.5-4.0 µm (e.g., F356W), please simply use the models for F322W2; for filters at 4.0-5.0µm (e.g., F410M), please use the models for F444W.

These models have a typical accuracy of RMS<0.1 pixel if astrometric and filter offsets are corrected. Below is a plot showing the $dy(dx=0)$ map for F444W module A grism R. The dy of my model is color-coded as the background image, and squares with filled colors is what I directly measured from the spectral traces in the LMC field. You can see a field variation of $dy(dx=0)$ of 1-2 pixel, which has been well modeled with my spectral tracing function.

dy(dx=0) map for F444W module A grism R

  1. Grism Dispersion model

The grism dispersion models give the relation of $dx(x_0, y_0, \lambda_s)$, i.e., the dispersion at given wavelength for source at $(x_0, y_0)$. In the commissioning, this was characterized by the F322W2 and F444W spectra of post-AGB star IRAS 05248-7007 in the LMC field.

The grism dispersion function is defined as fit_disp_order32() (or func_fit_wave(); see the beginning of this notebook). The best-fit parameters for these functions in different filter/module/grism can be found at http://u.arizona.edu/~fengwusun/data/dx_wave_dispersion_model/. The function for the same grism + module combination applies for all filters.

The accuracy of wavelength calibration is on the order of RMS=0.1-0.2 pixel. However, I would like to note that the current wavelength calibration data is still quite limited (~200 data points fit for a 16-parameter functions), and I would expect the accuracy gets further improved with cycle-1 calibration program.

  1. Sensitivity Function

All of the grism transmission curves that I have produced are available at http://u.arizona.edu/~fengwusun/data/sensitivity_model/. These files were derived using the Cycle-1 WFSS flux cal observations obtained in September, 2022. All of the grism R flux calibration is based on a G dwarf (P330-E; PID: 1538), a white dwarf (G191-B2B; PID: 1537) and an A-type star (J1743045; PID: 1536) with six exposures in total. The C grism observation of the white dwarf is corrupted by nearby stars, so most of the flux calibration of grism C is only based on P330E and J1743045 (4 individual integrations in total). The derived grism sensitivity functions in all filter+grism+module combinations are shown as follows:

overall sensitivity curve

The three columns of each grism sensitivity table are wavelength (in µm), sensitivity (in DN/s/Jy/pixel) and uncertainty of sensitivity. The sensitivity is the total throughput (OTE QE filter * grism, etc.) and sums the counts across the perpendicular (i.e. spatial) direction of the spectral trace, so it is in unit of DN/s/Jy/pixel. The dispersion is around 0.97 nm/pixel, so it is possible to convert the sensitivity to those in unit like (DN/s)/(erg/s/cm^2/Å)/nm. Note that the dispersion (nm/pixel or Å/pixel) is a function of wavelength and direct imaging position $(x0, y0)$, but 0.97 nm/pixel is a good characteristic dispersion.

Load the grism tracing/dispersion/sensitivity functions and POM-applied catalogs:

In [63]:
### select the filter that we will work on:
tmp_filter = 'F444W'
if tmp_filter == 'F444W': list_coords, tb_sex = coord_f444w, tb_f444w
if tmp_filter == 'F322W2': list_coords, tb_sex = coord_f322w2, tb_f322w2

### Interested wavelength Range
if tmp_filter == 'F444W': WRANGE = np.array([3.8, 5.1])
elif tmp_filter == 'F322W2': WRANGE = np.array([2.4, 4.1])
elif tmp_filter == 'F356W':  WRANGE = np.array([3.1, 4.0])
elif tmp_filter == 'F277W':  WRANGE = np.array([2.4, 3.1])

### Spectral tracing parameters:
if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
else: disp_filter = tmp_filter
tb_order23_fit_AR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'A', 'R'))
fit_opt_fit_AR, fit_err_fit_AR = tb_order23_fit_AR['col0'].data, tb_order23_fit_AR['col1'].data
tb_order23_fit_BR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'B', 'R'))
fit_opt_fit_BR, fit_err_fit_BR = tb_order23_fit_BR['col0'].data, tb_order23_fit_BR['col1'].data
tb_order23_fit_AC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'A', 'C'))
fit_opt_fit_AC, fit_err_fit_AC = tb_order23_fit_AC['col0'].data, tb_order23_fit_AC['col1'].data
tb_order23_fit_BC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'B', 'C'))
fit_opt_fit_BC, fit_err_fit_BC = tb_order23_fit_BC['col0'].data, tb_order23_fit_BC['col1'].data

### grism dispersion parameters:
tb_fit_displ_AR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('A', "R"))
w_opt_AR, w_err_AR = tb_fit_displ_AR['col0'].data, tb_fit_displ_AR['col1'].data
tb_fit_displ_BR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('B', "R"))
w_opt_BR, w_err_BR = tb_fit_displ_BR['col0'].data, tb_fit_displ_BR['col1'].data
tb_fit_displ_AC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('A', "C"))
w_opt_AC, w_err_AC = tb_fit_displ_AC['col0'].data, tb_fit_displ_AC['col1'].data
tb_fit_displ_BC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('B', "C"))
w_opt_BC, w_err_BC = tb_fit_displ_BC['col0'].data, tb_fit_displ_BC['col1'].data

### list of module/pupil and corresponding tracing/dispersion function:
list_mod_pupil   = np.array(['AR', 'BR', 'AC', 'BC'])
list_fit_opt_fit = np.array([fit_opt_fit_AR, fit_opt_fit_BR, fit_opt_fit_AC, fit_opt_fit_BC])
list_w_opt       = np.array([w_opt_AR, w_opt_BR, w_opt_AC, w_opt_BC])

### Sensitivity curve:
dir_fluxcal = '/home/u24/fengwusun/Comissioning/cycle1_cal_program/fluxcal/product/'
tb_sens_AR = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat' % (disp_filter, 'A', 'R'))
tb_sens_BR = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat'% (disp_filter, 'B', 'R'))
tb_sens_AC = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat' % (disp_filter, 'A', 'C'))
tb_sens_BC = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat'% (disp_filter, 'B', 'C'))
f_sens_AR = interpolate.UnivariateSpline(tb_sens_AR['wavelength'], tb_sens_AR['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_BR = interpolate.UnivariateSpline(tb_sens_BR['wavelength'], tb_sens_BR['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_AC = interpolate.UnivariateSpline(tb_sens_AC['wavelength'], tb_sens_AC['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_BC = interpolate.UnivariateSpline(tb_sens_BC['wavelength'], tb_sens_BC['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
list_f_sens       = (f_sens_AR, f_sens_BR, f_sens_AC, f_sens_BC)

'''list of source list (POM-applied catalog)'''
list_sl = glob.glob('/home/u24/fengwusun/IceAge/POM_catalog/jw01309*_nrc[a-b]_%s_dirimg_sources.list' % tmp_filter)
list_sl.sort()
list_rate_this = np.array([os.path.dirname(all_rate_list[0]) + '/' + os.path.basename(list_sl[i]).split('_%s' % tmp_filter)[0] + 'long_rate_lv1.5.fits' 
                           for i in range(len(list_sl))])

### apply magnitude cut to the input catalog?
mag_cut = 18
Index_interest = tb_sex[tb_sex['abmag'] <= mag_cut]['id'].data
# Index_interest = tb_sex[tb_sex['MAG_AUTO'] <= mag_cut]['NUMBER'].data
# Index_interest = tb_sex[tb_sex['FLUX_AUTO'] >= 10**(-0.4 * (mag_cut - nircam_LW_ZP))]['NUMBER'].data
print(len(Index_interest), 'sources selected')
list_sl_table = [ascii.read(x) for x in list_sl]


''' 
Generate two dictionaries which record the grism frame ID and path for each source to be extracted
'''
grism_frame_ID_per_source = dict()   # dictionary of the ID of exposure (0-95) for each source (0-5000+)
grism_frame_path_per_source = dict() # dictionary of the path of exposure (*/jw*.fits) for each source (0-5000+)
for idx in Index_interest[:]:
    grism_frame_ID_per_source['%s' % idx] = []
    grism_frame_path_per_source['path_%s' % idx] = []
for i, tmp_sl_table in enumerate(list_sl_table):
    for idx in tmp_sl_table['Index'].data:
        try:
            grism_frame_ID_per_source['%s' % idx].append(i)
            grism_frame_path_per_source['path_%s' % idx].append(list_rate_this[i])
        except KeyError: continue
N_frame_per_source = np.array([len(grism_frame_ID_per_source[x]) for x in grism_frame_ID_per_source.keys()])
9 sources selected

The following functions are used to extract 2d spectra, resample / coadd 2d spectra and save them.

In [66]:
def extract_2d_spec(img, WRANGE, x0, y0, dxs, dys, wave, img_wht, img_dq, header, 
                    aper = 10., pupil = 'R', n_procs = 4):
    '''
    Extract 2D spectrum from the Grism Image
    -----------------------------------------------
        Parameters
        ----------
        img : `~numpy.ndarray`
            2D dispersed slitless spectroscopic image
        
        GConf: `~grismconf.grismconf.Config`
            Configuration file of grism (taken from `grismconf`)
            
        x0, y0 : float
            Reference position (i.e., in direct image)
        
        dxs, dys : `~numpy.ndarray`
            First derivative of the trace x/y-coordinates with respect to 't' (taken from `grismconf`)
        
        wave: 
            Wavelength of extracted spectrum
            
        img_wht : `~numpy.ndarray`
            Weighting image (1/rms^2) of `img`
            
        img_dq : `~numpy.ndarray`
            DQ (data-quality) image (1/rms^2) of `img`
            
        header: `fits.Header()` object
            Header of 2D dispersed slitless spectroscopic image
            
        aper: float
            Aperture radius for 2D spectra extraction (unit: arcsec, default: 1.25)
        
        pix_scale: float
            Pixel scale in unit of arcsec. If None then computed from header
        
        pupil: 'R' or 'C'
            pupil of grism ('R' or 'C')
            
        Returns
        -------
        hdul : `~astropy.io.fits.HDUList()` object
            The fits files of extracted 2D spectra, including a primary frame, a sci frame,
                (scientific), a wht frame (weighting) and a dq frame (data quality).
    '''
    w_min, w_max = WRANGE
    x_on_G_img = dxs + x0
    y_on_G_img = dys + y0
    
    if pupil not in ['C', 'R']: raise KeyError('pupil of grism should be either "R" or "C"! ')
    
    aper_int_pix = int(aper)
    
    
    ### region that we can extract spectra: 
    if pupil == 'R':
        args_eff = np.where((wave >= w_min) & (wave <= w_max) 
                            & (x_on_G_img >= 5) & (x_on_G_img <= 2047 - 6)
                            & (y_on_G_img >= aper_int_pix) & (y_on_G_img <= 2047 - 6 - aper_int_pix) )
    elif pupil == 'C':
        args_eff = np.where((wave >= w_min) & (wave <= w_max) 
                            & (x_on_G_img >= aper_int_pix) & (x_on_G_img <= 2047 - 6 - aper_int_pix)
                            & (y_on_G_img >= 5) & (y_on_G_img <= 2047 - 6) )
    if np.size(args_eff) <= 20:  ## spectrum is too short
        raise ValueError('No spectrum can be extracted from the 2D Image')
    tmp_spec_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
    tmp_wht_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
    tmp_dq_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
    
    for i, j in enumerate(args_eff[0]):
        if pupil == 'R':
            tmp_x, tmp_y1, tmp_y2 = int(x_on_G_img[j]), int(y_on_G_img[j] - aper_int_pix - 1), int(y_on_G_img[j] + aper_int_pix + 2)
            img.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
            img_wht.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img_wht.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
            img_dq.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img_dq.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
            tmp_spec_2d[i] = img.T[tmp_x][tmp_y1+1:tmp_y2-1]
            tmp_wht_2d[i] = img_wht.T[tmp_x][tmp_y1+1:tmp_y2-1]
            tmp_dq_2d[i] = img_dq.T[tmp_x][tmp_y1+1:tmp_y2-1]
        elif pupil == 'C':
            tmp_y, tmp_x1, tmp_x2 = int(y_on_G_img[j]), int(x_on_G_img[j] - aper_int_pix - 1), int(x_on_G_img[j] + aper_int_pix + 2)
            img[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
            img_wht[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img_wht[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
            img_dq[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img_dq[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
            tmp_spec_2d[i] = img[tmp_y, tmp_x1+1:tmp_x2-1]
            tmp_wht_2d[i] = img_wht[tmp_y, tmp_x1+1:tmp_x2-1]
            tmp_dq_2d[i] = img_dq[tmp_y, tmp_x1+1:tmp_x2-1]
    
    tmp_spec_2d = tmp_spec_2d.T
    tmp_wht_2d = tmp_wht_2d.T
    tmp_dq_2d = tmp_dq_2d.T
        
    '''Construct fits file for 2d grism spectra'''
    ### Primary HDU
    hdu = fits.PrimaryHDU()
    hdu.header['x0'] = (np.float32(x0), 'Reference position X in direct image')
    hdu.header['y0'] = (np.float32(y0), 'Reference position Y in direct image')
    # tmp_coord = wcs.utils.pixel_to_skycoord(x0, y0, tmp_wcs)
    # hdu.header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
    # hdu.header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
    hdu.header['author'] = ('Fengwu Sun', 'Author of this file')
    hdu.header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S",  time.localtime()), 'Time of Creation')
    ### Scientific HDU
    hdu_sci = fits.ImageHDU(np.float32(tmp_spec_2d), name = 'SPEC2D')
    hdu_sci.header['wave_1'] = (wave[args_eff[0][1]], 'Wavelength (um) of first pixel')
    hdu_sci.header['d_wave'] = (np.mean(np.diff(wave[args_eff[0]])), 'Wavelength Difference (um) between each pixel')
    hdu_sci.header['comments'] = ('wave = wave_1 + np.arange(0, NAXIS1, 1) * d_wave')
    # hdu_sci.header['pixscale'] = (pix_scale, 'Pixel scale in undispersed direction (arcsec)')
    hdu_sci.header['aperture'] = (aper, 'Aperture radius in undispersed direction (arcsec)')
    hdu_sci.header['pupil'] = (pupil, 'Pupil of grism (R or C)')
    hdu_sci.header['module'] = (header['module'], 'Module of Detector (A or B)')
    hdu_sci.header['diff_y'] = (y_on_G_img[args_eff[0]][0], 'Y_(full)-Y_(trim)')
    hdu_sci.header['diff_x'] = (x_on_G_img[args_eff[0]][0], 'X_(full)-X_(trim)')
    ### Weight / Data-quality HDU
    hdu_wht = fits.ImageHDU(np.float32(tmp_wht_2d), name = 'WHT2D')
    hdu_dq = fits.ImageHDU(np.int32(tmp_dq_2d), name = 'DQ2D')
    ### Tracing & Dispersion information
    tmp_tb_wave = Table(data = [wavs[args_eff], x_on_G_img[args_eff], y_on_G_img[args_eff], dxs[args_eff], dys[args_eff]],
                    names = ['wavelength', 'xs', 'ys', 'dxs', 'dys'])
    tmp_tb_wave['wavelength'].info.format = '.6f'
    for x in tmp_tb_wave.colnames[1:]: tmp_tb_wave[x].info.format = '.3f'
    hdul = fits.HDUList([hdu, hdu_sci, hdu_wht, hdu_dq])
    hdul.append(fits.BinTableHDU(tmp_tb_wave, name = 'WAVE'))
    return hdul


def store_all_2d_spec(fits_list, pupils, modules, paths, output = 'test.fits',
                      coord = None, filter = None, info_table = None, overwrite = True):
    '''
    Write all the individual 2D spectra of one single sources into a fits file
    -----------------------------------------------
    Parameters
        ----------  
        fits_list : list of `astropy.io.fits`
            list of fits files generated by extract_2d_spec()
        
        pupils: list of strings
            pupil list of grism ('R' or 'C')
        
        modules: list of strings
            module list of NIRCam detector ('A' or 'B')
        
        paths: list of strings
            paths of grism exposures
        
        # poms: list of integers
        #     list of POM values
        
        output: string
            Output file name
        
        coord: `astropy.coordinates.SkyCoord` object
            coordinates of the source. default: None (will not record this in the header)
        
        filter: string
            Filter of grism observation (e.g., 'F444W', 'F332W', 'F356W')
            default: None (will not record this in the header)
        
        info_table: astropy.table.Row
            A row of information related to this source
            default: None (will not record this in the header)
            
        overwrite: bool
            If true, save (and overwrite) the data. default: True
            
        Returns
        -------
        ind_hdul: `~astropy.io.fits.HDUList`
            fits HDU list of individual sci, wht and dq data.
            
    '''
    ## HDU list of individual exposures:
    ind_hdul = fits.HDUList([fits.PrimaryHDU()])
    for l, x in enumerate(fits_list):
        ind_hdul.append(x[1]) 
        ind_hdul[-1].header['EXTNAME'] = 'SPEC2D-%d' % l
        for card_name in ['x0', 'y0']: # 'ra0', 'dec0']:
            ind_hdul[-1].header[card_name] = x[0].header[card_name]
        ind_hdul[-1].header['pupil'] = pupils[l]
        ind_hdul[-1].header['module'] = modules[l]
        # ind_hdul[-1].header['pom'] = poms[l]
        ind_hdul[-1].header['datapath'] = paths[l].split('/')[-1]
        ind_hdul.append(x[2]) 
        ind_hdul[-1].header['EXTNAME'] = 'WHT2D-%d' % l
        ind_hdul.append(x[3])
        ind_hdul[-1].header['EXTNAME'] = 'DQ2D-%d' % l
        ind_hdul.append(x[4])
        ind_hdul[-1].header['EXTNAME'] = 'WAVE-%d' % l
    #### stats of single exposures
    ind_hdul.append(fits.BinTableHDU(Table(names = ['id', 'pupil', 'module', 'datapath'], # 'POM', 
                                           data = [range(len(fits_list)), pupils, modules, # poms,
                                                   [tmp_path.split('/')[-1] for tmp_path in paths]]),
                                     name = 'STATS'))
    #### add other important specs
    ind_hdul[0].header['DIRNAME'] = (os.path.dirname(paths[0]), 'Directory name of original grism data')
    if coord != None:
        ind_hdul[0].header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
        ind_hdul[0].header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
    ind_hdul[0].header['N_coadd'] = (len(ind_hdul[-1].data['pupil']), 'Number of coadded frames')
    ind_hdul[0].header['N_R'] = (np.sum(ind_hdul[-1].data['pupil'] == 'R'), 'Number of R grism frames')
    ind_hdul[0].header['N_C'] = (np.sum(ind_hdul[-1].data['pupil'] == 'C'), 'Number of C grism frames')
    ind_hdul[0].header['author'] = ('Fengwu Sun', 'Author of this file')
    ind_hdul[0].header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S",  time.localtime()), 'Time of Creation')
    if filter != None:
        ind_hdul[0].header['filter'] = (tmp_filter, 'Filter name')
    if info_table != None:
        ind_hdul[0].header['COMMENTS'] = 'Belows are information taken from input catalog:'
        for x in info_table.colnames:
            ind_hdul[0].header['HIERARCH ' + x] = info_table[x]
    if overwrite: ind_hdul.writeto(output, overwrite = overwrite) ## save the fits file?
    return ind_hdul

def resample_spec2d_wmin_wmax(x, i):
    '''
    Resample a series of 2D spectra at wavelegth range = tmp_w_min - tmp_w_max
    -----------------------------------------------
    Parameters
        ---------- 
        x: int
            index of the individual fits files in the list.
        
        i: int
            index of the resampled wave array
        
        tmp_wave_2d: 
            Global variable (`list(np.array())`). 
            A compilation of 1D wavelength array associated with all individual grism frames.
        
        tmp_ind_fits_list:
            Global variable (`astropy.io.fits.HDUList()`). 
            A HDU List that have extracted 2D spectra (SCI/WHT/DQ) from each individual frames
            
        wave_sample:
            Global variable (`np.array()`)
            Resampled wavelengths. 
            
        Returns
        -------
        tmp_spec_w_1: `~np.array()`
            a column of 2D SCI spectra from frame [x] that fall in the resampled wavelength range [i]
        
        tmp_wht_w_1: `~np.array()`
            a column of 2D WHT spectra at the same position as that of `tmp_spec_w_1`
        
        tmp_cov_w_1: `~np.array()`
            a column of 2D DQ spectra at the same position as that of `tmp_spec_w_1`
    '''
    global tmp_wave_2d, tmp_ind_fits_list, wave_sample # arrays to use
    # global arr_spec_2d, arr_wht_2d, arr_cov_2d         # arrays to save; impossible with pool
    tmp_w_min, tmp_w_max = wave_sample[i], wave_sample[i+1]
    ## resampled SCI, WHT, DQ image
    arg_in_wmin_wmax = tuple([(tmp_wave_2d[x] > tmp_w_min) & (tmp_wave_2d[x] <= tmp_w_max)])
    tmp_spec_w = tmp_ind_fits_list[x][1].data.T[arg_in_wmin_wmax]
    tmp_wht_w = np.nan_to_num(tmp_ind_fits_list[x][2].data.T[arg_in_wmin_wmax], posinf = 0, neginf = 0)
    tmp_dq_w = tmp_ind_fits_list[x][3].data.T[arg_in_wmin_wmax]
    tmp_wht_w[tmp_dq_w!=0] = 0
    ## stack/resample N spectra at wmin-wmax (M rows)
    tmp_spec_w_1 = np.nansum(tmp_spec_w * tmp_wht_w, axis = 0) / np.nansum(tmp_wht_w, axis = 0) #/ sum_tmp_wht_w #
    tmp_wht_w_1 = np.nansum(tmp_wht_w, axis = 0)
    tmp_cov_w_1 = np.int8(tmp_wht_w_1 != 0)
    ## save resampled spectra to array
    # arr_spec_2d[i, x], arr_wht_2d[i, x], arr_cov_2d[i, x] = tmp_spec_w_1, tmp_wht_w_1, tmp_cov_w_1
    return(tmp_spec_w_1, tmp_wht_w_1, tmp_cov_w_1)
In [ ]:
 

Extract 2D grism spectra with the loop below:

In [68]:
'''Choose extraction mode:
all : R + C combined
  R : R only
  C : C only
'''

extract_mode = 'all'   # all, R, C
start = time.time()

for i, idx in enumerate(Index_interest[:]):
    '''
    Read source information from the input catalog 
    '''
    t_a = time.time()
    ## item_sex = tb_sex[tb_sex['NUMBER'] == idx][0]
    ## tmp_mag = item_sex['MAG_AUTO']
    item_sex = tb_sex[tb_sex['id'] == idx][0]
    tmp_mag = item_sex['abmag']
    tmp_coord = SkyCoord(item_sex['RA'], item_sex['DEC'], unit = (u.deg, u.deg))
    ### tmp_coord = SkyCoord(item_sex['X_WORLD'], item_sex['Y_WORLD'], unit = (u.deg, u.deg))
    print('-[%3d] ID%-5d mag=%.2f' % (i, idx, tmp_mag), end = " ")
    '''
    ### if no record of grism exposure is found in the dict `grism_frame_ID_per_source`, 
    ###    the script will continue to save time
    '''
    if len(grism_frame_ID_per_source[str(idx)]) == 0: print(' >> no spec found; '); continue
    else:
        ## find 2D grism image, get the list
        this_list_rate_wcs = np.array(grism_frame_path_per_source['path_'+str(idx)])
        this_ID_rate_wcs = np.array(grism_frame_ID_per_source[str(idx)])  
    # print('[%d spectra]' % len(this_list_rate_wcs))
    ### but, if the number of available frame is too small, we will skip the extraction
    if len(grism_frame_ID_per_source[str(idx)]) <= 2: print(' >> no spec found; '); continue
    
    
    '''
    Examine each reduced image to extract 2D spectra
    '''
    tmp_ind_fits_list = []
    tmp_ind_fits_name_list = []
    tmp_ind_module_list = []
    tmp_ind_pupil_list = []
    tmp_ind_pom_list = []
    for j, tmp_grism in enumerate(this_list_rate_wcs):
        ### `tmp_grism` is the path of grism data ('*/jw*.fits')
        # if idx not in list_sl_table[j]['Index'] : continue  ## legacy command
        ### `tmp_idx_grism` is the sequential number of grism data (like 0-95 with 96 frames)
        tmp_idx_grism = this_ID_rate_wcs[j]  
        tmp_grism_img = fits.getdata(tmp_grism, 'sci')  ### Grism SCI image data
        tmp_grism_dq = fits.getdata(tmp_grism, 'dq')    ### Grism DQ data 
        tmp_grism_hd = fits.getheader(tmp_grism)        ### Grism primary Header
        tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
        ### depending on `extract_mode`, only extrac R/C data
        if extract_mode.lower() == 'r': 
            if tmp_pupil == 'C' : continue
        if extract_mode.lower() == 'c': 
            if tmp_pupil == 'R' : continue
        tmp_grism_hd_sci = fits.getheader(tmp_grism, 'sci') ### Grism SCI Header
        ### Read or Generate Grism WHT data:
        tmp_grism_wht_path = tmp_grism.replace('lv1.5.fits', 'wht.fits')
        if os.path.isfile(tmp_grism_wht_path): 
            tmp_grism_wht = fits.getdata(tmp_grism_wht_path)
        else:
            ### 2D grism WHT file not found, generate from Error map:
            tmp_grism_err = fits.getdata(tmp_grism, 'err')
            tmp_grism_err[tmp_grism_err == 0] = np.nan
            tmp_grism_wht = tmp_grism_err**-2
            tmp_grism_wht_hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(tmp_grism_wht, name = 'WHT')])
            tmp_grism_wht_hdul[1].header = tmp_grism_hd_sci
            tmp_grism_wht_hdul.writeto(tmp_grism_wht_path, overwrite = True)
            tmp_grism_wht_hdul.close()
        ### Direct imaging positions 
        item_sl = list_sl_table[tmp_idx_grism][list_sl_table[tmp_idx_grism]['Index'] == idx][0]
        x0 = item_sl['pixel_x']
        y0 = item_sl['pixel_y']
        print('%s(%.1f, %.1f) ' % (tmp_pupil, x0, y0), end = ' ')
        # if (x0 > 2100) & (idx == 1357) & (tmp_pupil == 'R'): continue
        ### spectral tracing parameters
        fit_opt_fit = list_fit_opt_fit[list_mod_pupil == tmp_module + tmp_pupil][0]
        w_opt = list_w_opt[list_mod_pupil == tmp_module + tmp_pupil][0]
        dxs, dys, wavs = grism_conf_preparation(x0 = x0, y0 = y0, pupil = tmp_pupil, 
                                                fit_opt_fit = fit_opt_fit, w_opt = w_opt)
        tmp_aper = 15.  # aperture size 
        ### Extract 2D spectrum on each single frame
        try:
            tmp_ind_fits = extract_2d_spec(img = tmp_grism_img, WRANGE = WRANGE, 
                                           x0 = x0, y0 = y0, dxs = dxs, dys = dys,
                                           wave = wavs, img_wht = tmp_grism_wht, img_dq = tmp_grism_dq,
                                           header = tmp_grism_hd, pupil = tmp_pupil, aper = tmp_aper)
        except ValueError:
            continue
        tmp_ind_fits_list.append(tmp_ind_fits)
        tmp_ind_module_list.append(tmp_module)
        tmp_ind_fits_name_list.append(tmp_grism)
        tmp_ind_pupil_list.append(tmp_pupil)
        
    ### If no spectra can be extract, continue; else: save spectra of all individal extraction
    if len(tmp_ind_fits_list) == 0: 
        print(' >> no spec found; ')
        continue
    else: 
        print(' >> N=%d, x0=%.1f, y0=%.1f' % (len(tmp_ind_fits_list), x0, y0), end = '')
        if extract_mode.lower() not in ['r', 'c']:
            '''Change the path of saved 2D spectra: '''
            tmp_ind_specs_name = '/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/allspec_2d_%s_%s_ID%s.fits' % (tmp_target, tmp_filter, idx)
            store_all_2d_spec(fits_list = tmp_ind_fits_list, pupils = tmp_ind_pupil_list, 
                              modules = tmp_ind_module_list, paths = tmp_ind_fits_name_list, 
                              output = tmp_ind_specs_name, coord = tmp_coord, 
                              filter = tmp_filter,  info_table = item_sex, overwrite = True)
            print(' >> save all extracted spec2d ', end = '')
    ### Correct for sensitivity:
    for k, tmp_fits in enumerate(tmp_ind_fits_list[:]):
        tmp_modpup = tmp_fits[1].header['MODULE'] + tmp_fits[1].header['PUPIL']
        tmp_f_sens = list_f_sens[np.where(list_mod_pupil == tmp_modpup)[0][0]]
        tmp_wavelength = tmp_fits[4].data['wavelength']
        tmp_ind_fits_list[k][1].header['bunit'] = ('mJy', 'Brightness Unit')
        tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy
        tmp_ind_fits_list[k][2].data = tmp_fits[2].data * tmp_f_sens(tmp_wavelength)**2 * 1e-6
    '''
    Combine 2D spectra from all extraction
    '''
    tmp_wave_2d = np.array([x[-1].data['wavelength'] for x in tmp_ind_fits_list], dtype = object)
    ### resampled to the following wavelength frame:
    d_wave = 0.001
    wave_sample = np.arange(WRANGE[0], WRANGE[1] + d_wave, d_wave)
    wave_sample_c = wave_sample[:-1] + d_wave / 2.
    arr_spec_2d = np.zeros((len(wave_sample)-1, len(tmp_wave_2d), len(tmp_ind_fits_list[0][1].data)))
    arr_wht_2d = np.zeros_like(arr_spec_2d)
    arr_cov_2d = np.zeros_like(arr_spec_2d)
    np.seterr(invalid = 'ignore')
    ### New method: resample spectrum using a multiple-processing tool:
    xx, ii = np.meshgrid(range(len(tmp_wave_2d)), range(len(wave_sample) - 1))
    xx = xx.flatten()
    ii = ii.flatten()
    with Pool(np.min([8, len(xx)])) as pool:
        output = pool.starmap(resample_spec2d_wmin_wmax,  zip(xx, ii))
    for o in range(len(output)):
        i, x = ii[o], xx[o] 
        arr_spec_2d[i, x], arr_wht_2d[i, x], arr_cov_2d[i, x] = output[o]
    ## sigma clip spectra SED 
    tmp_sigma_clip = SigmaClip(sigma_lower = 2.0, sigma_upper = 3.0)
    arg_sigclip = np.where(tmp_sigma_clip(np.nanmean(arr_spec_2d, axis = (0, 2))).mask == False)[0]
    tmp_ind_fits_list = [tmp_ind_fits_list[x] for x in arg_sigclip]
    tmp_ind_module_list = np.array(tmp_ind_module_list)[arg_sigclip]
    tmp_ind_fits_name_list = np.array(tmp_ind_fits_name_list)[arg_sigclip]
    arr_spec_2d = arr_spec_2d[:,arg_sigclip,:]
    arr_wht_2d = arr_wht_2d[:,arg_sigclip,:]
    arr_cov_2d = arr_cov_2d[:,arg_sigclip,:]
    ## weighted mean 2d spectra:
    tmp_spec_2d = np.nansum(arr_spec_2d * arr_wht_2d, axis = 1) / np.nansum(arr_wht_2d, axis = 1) # sum_arr_wht_2d # 
    tmp_wht_2d = np.nansum(arr_wht_2d, axis = 1)
    tmp_cov_2d = np.nansum(arr_cov_2d, axis = 1)
    '''
    Save Coadded 2D spectra
    '''
    tmp_tb_cov = Table(names = ['index', 'name', 'x0', 'y0', 'module', 'pupil', 'DIFF_X', 'DIFF_Y', 'wave_0', 'wave_1'],
                       data = [range(len(tmp_ind_fits_list)),
                               [tmp_grism.split('/')[-1][:-16] for tmp_grism in tmp_ind_fits_name_list],
                               [x[0].header['x0'] for x in tmp_ind_fits_list],
                               [x[0].header['y0'] for x in tmp_ind_fits_list],
                               tmp_ind_module_list,
                               [x[1].header['PUPIL'] for x in tmp_ind_fits_list],
                               [x[1].header['DIFF_X'] for x in tmp_ind_fits_list],
                               [x[1].header['DIFF_Y'] for x in tmp_ind_fits_list],
                               [np.round(wave_sample_c[np.sum(arr_wht_2d, axis = -1)[:, n_] != 0][0], 4) 
                                for n_ in range(len(tmp_ind_fits_list))],
                               [np.round(wave_sample_c[np.sum(arr_wht_2d, axis = -1)[:, n_] != 0][-1], 4) 
                                for n_ in range(len(tmp_ind_fits_list))],
                              ]) 
    
    hdu = fits.PrimaryHDU()
    hdu.header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
    hdu.header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
    hdu.header['N_coadd'] = (len(tmp_tb_cov['pupil']), 'Number of coadded frames')
    hdu.header['N_R'] = (np.sum(tmp_tb_cov['pupil'] == 'R'), 'Number of R grism frames')
    hdu.header['N_C'] = (np.sum(tmp_tb_cov['pupil'] == 'C'), 'Number of C grism frames')
    hdu.header['author'] = ('Fengwu Sun', 'Author of this file')
    hdu.header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S",  time.localtime()), 'Time of Creation')
    hdu.header['filter'] = (tmp_filter, 'Filter name')
    # hdu.header['abmag'] = (tmp_mag, 'AB Magnitude in %s Band' % tmp_filter)
    hdu.header['COMMENTS'] = 'Belows are information taken from input mock catalog:'
    for x in item_sex.colnames:
        hdu.header['HIERARCH ' + x] = item_sex[x]
    hdu_sci = fits.ImageHDU(tmp_spec_2d.T, name = 'SPEC2D')
    hdu_sci.header['wave_1'] = (wave_sample_c[0], 'Wavelength (um) of first pixel')
    hdu_sci.header['d_wave'] = (wave_sample_c[1] - wave_sample_c[0], 'Wavelength Difference (um) between each pixel')
    hdu_sci.header['comments'] = ('wave = wave_1 + np.arange(0, NAXIS1, 1) * d_wave')
    hdu_sci.header['pixscale'] = (0.0629, 'Pixel scale in undispersed direction (arcsec)')
    hdu_sci.header['aperture'] = (tmp_aper, 'Aperture radius in undispersed direction (arcsec)')
    hdu_wht = fits.ImageHDU(tmp_wht_2d.T, name = 'WHT2D')
    hdu_wht.header['comments'] = ('Weight image; ERR = WHT^(-0.5)')
    hdu_cov = fits.ImageHDU(tmp_cov_2d.T, name = 'COV2D')
    hdu_cov.header['comments'] = ('Coverage image')
    hdu_tab = fits.BinTableHDU(tmp_tb_cov, name = 'STATS')
    hdu_tab.header['COMMENT'] = 'name:     name of simulated image'
    hdu_tab.header['COMMENT'] = 'x0 / y0:  reference position (i.e., in direct image)'
    hdu_tab.header['COMMENT'] = '   (in reality this is registered to the wcs of grism image)'
    hdu_tab.header['COMMENT'] = 'DIFF_X:   X_(full)-X_(trim)'
    hdu_tab.header['COMMENT'] = 'DIFF_Y:   Y_(full)-Y_(trim)'
    hdu_tab.header['COMMENT'] = 'wave_0:   minimum wavelength (micron) in this coverage'
    hdu_tab.header['COMMENT'] = 'wave_1:   maximum wavelength (micron) in this coverage'
    hdul = fits.HDUList([hdu, hdu_sci, hdu_wht, hdu_cov, hdu_tab])
    '''Change the path of saved 2D spectra: '''
    tmp_spec_2d_name = '/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/spec_2d_%s_%s_ID%s_comb.fits' % (tmp_target, tmp_filter,  idx)
    if extract_mode.lower() in ['r', 'c']: 
        tmp_spec_2d_name = tmp_spec_2d_name.replace('_comb.fits', '_%s.fits' % extract_mode.upper())
    hdul.writeto(tmp_spec_2d_name, overwrite = True)
    print(' >> save stacked spec2d; ')

end = time.time()
print('run time = %.1fs = %.1fmin' % (end - start, (end - start)/60))
-[  0] ID521   mag=16.70 R(15.5, 1304.2)  R(25.1, 1302.4)  R(13.3, 1313.8)  R(6.1, 1294.6)  R(136.4, 1186.3)  R(146.0, 1184.3)  R(134.4, 1195.7)  R(126.9, 1176.7)  R(257.2, 1068.2)  R(266.8, 1066.3)  R(255.1, 1077.7)  R(247.7, 1058.6)   >> N=12, x0=247.7, y0=1058.6
2022-10-30 23:04:39,360 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
-[  1] ID498   mag=16.86  >> no spec found; 
-[  2] ID805   mag=16.87 R(1087.0, 1484.9)  R(1096.6, 1482.9)  R(1085.0, 1494.5)  R(1077.5, 1475.4)  R(1206.9, 1365.5)  R(1216.4, 1363.5)  R(1204.9, 1374.9)  R(1197.3, 1356.0)  R(1326.3, 1246.1)  R(1335.9, 1244.1)  R(1324.3, 1255.5)  R(1316.8, 1236.5)  R(1767.8, 1493.9)  R(1777.2, 1491.8)  R(1765.8, 1503.4)  R(1758.1, 1484.5)  R(1886.0, 1373.9)  R(1895.5, 1371.7)  R(1884.1, 1383.3)  R(1876.4, 1364.5)  R(2003.8, 1254.2)  R(2013.3, 1252.2)  R(2001.9, 1263.6)  R(1994.2, 1244.9)  C(1109.8, 1089.5)  C(1119.3, 1087.6)  C(1107.8, 1098.9)  C(1100.2, 1080.0)  C(1229.2, 970.7)  C(1238.7, 968.7)  C(1227.2, 980.2)  C(1219.7, 961.3)  C(1348.3, 852.2)  C(1357.7, 850.2)  C(1346.2, 861.6)  C(1338.8, 842.8)  C(1109.8, 1089.5)  C(1119.3, 1087.6)  C(1107.8, 1098.9)  C(1100.2, 1080.0)  C(1229.2, 970.7)  C(1238.7, 968.7)  C(1227.2, 980.1)  C(1219.7, 961.2)  C(1348.3, 852.2)  C(1357.7, 850.2)  C(1346.2, 861.6)  C(1338.8, 842.8)   >> N=48, x0=1338.8, y0=842.8
2022-10-30 23:04:49,601 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
-[  3] ID1540  mag=17.23 R(355.4, 69.1)  R(364.8, 67.2)  R(353.4, 78.4)  
2022-10-30 23:04:52,471 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

R(346.1, 59.9)   >> N=4, x0=346.1, y0=59.9 >> save all extracted spec2d  >> save stacked spec2d; 
-[  4] ID563   mag=17.26 C(284.0, 1917.4)  C(293.5, 1915.6)  C(282.0, 1927.1)  C(274.5, 1907.9)  C(284.0, 1917.4)  C(293.6, 1915.6)  C(281.9, 1927.0)  
2022-10-30 23:04:54,541 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

C(274.5, 1907.9)   >> N=8, x0=274.5, y0=1907.9 >> save all extracted spec2d  >> save stacked spec2d; 
-[  5] ID917   mag=17.34 R(48.0, 1745.6)  R(57.5, 1743.7)  R(45.9, 1755.1)  R(38.6, 1736.0)  R(168.5, 1628.0)  R(178.1, 1626.1)  R(166.4, 1637.5)  R(159.1, 1618.5)  R(289.0, 1510.2)  R(298.5, 1508.3)  R(286.8, 1519.7)  R(279.5, 1500.7)   >> N=12, x0=279.5, y0=1500.7
2022-10-30 23:04:57,248 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
-[  6] ID1509  mag=17.34 R(1560.0, 961.0)  R(1569.5, 959.0)  R(1558.1, 970.6)  R(1550.4, 951.6)  R(1678.4, 841.4)  R(1687.8, 839.4)  R(1676.4, 850.9)  R(1668.7, 832.1)  R(1796.1, 721.8)  R(1805.6, 719.7)  R(1794.2, 731.3)  R(1786.6, 712.4)  C(1554.2, 428.3)  C(1563.6, 426.3)  C(1552.3, 437.7)  C(1544.7, 419.0)  C(1671.3, 309.7)  C(1680.8, 307.6)  C(1669.4, 319.0)  C(1661.7, 300.3)  C(1788.2, 191.4)  C(1797.5, 189.3)  C(1786.1, 200.6)  C(1778.5, 182.1)  C(1554.2, 428.3)  C(1563.6, 426.2)  C(1552.4, 437.7)  C(1544.7, 419.0)  C(1671.4, 309.7)  C(1680.7, 307.6)  C(1669.4, 319.0)  C(1661.8, 300.4)  C(1788.2, 191.4)  C(1797.5, 189.3)  C(1786.1, 200.6)  C(1778.5, 182.0)   >> N=36, x0=1778.5, y0=182.0
2022-10-30 23:05:05,319 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
-[  7] ID1619  mag=17.74 R(1112.2, 214.7)  R(1121.6, 212.7)  R(1110.3, 224.1)  R(1102.8, 205.5)  R(1229.7, 97.1)  R(1239.0, 95.1)  R(1227.7, 106.4)  R(1220.2, 87.9)  R(1783.3, 221.0)  R(1792.6, 218.9)  R(1781.4, 230.4)  R(1773.7, 211.8)  R(1899.3, 102.6)  R(1908.8, 100.6)  R(1897.6, 112.0)  R(1889.8, 93.4)   >> N=16, x0=1889.8, y0=93.4
2022-10-30 23:05:09,288 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
-[  8] ID1351  mag=17.79 R(1275.8, 200.9)  R(1285.2, 199.0)  R(1273.8, 210.3)  R(1266.4, 191.6)  R(1393.7, 84.5)  R(1403.0, 82.5)  R(1391.7, 93.7)  R(1384.3, 75.1)  R(1944.8, 213.1)  R(1954.1, 211.1)  R(1942.9, 222.5)  R(1935.3, 203.8)   >> N=12, x0=1935.3, y0=203.8
2022-10-30 23:05:11,579 - stpipe - WARNING - /tmp/ipykernel_16422/919471653.py:107: RuntimeWarning: divide by zero encountered in true_divide
  tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy

 >> save all extracted spec2d  >> save stacked spec2d; 
run time = 35.6s = 0.6min
In [ ]:
 
In [ ]:
 

7. Examine 2D & 1D spectra in plots

With the 2D grism spectra you just extracted and combined above, it is now the time to examine the quality of extracted spectra. You can plot the 2D, 1D spectrum and direct imaging data in the same frame.

In [78]:
'''List of 2D co-added grism spectra that we just extracted:'''
list_spec = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/spec_2d_CHAMMS1-CO_F444W_ID*_comb.fits'))
print(len(list_spec), 'spectra files in total.')

'''Direct imaging data 
   for IceAge data, it is avaialble at:
       http://gxn.as.arizona.edu/~sunfengwu/IceAge/mosaics/IceAge_CHAMMS1-C2-FIELD_lw_F430M_visitall_modall_i2d.fits.tar.gz
'''
path_LW = '/xdisk/egami/fengwusun/ERSs/IceAge/all_mosaics/IceAge_CHAMMS1-C2-FIELD_lw_F430M_visitall_modall_i2d.fits'
img_LW = fits.getdata(path_LW, 'sci')
hd_LW =  fits.getheader(path_LW, 'sci')
wcs_LW = wcs.WCS(hd_LW)
pix_LW = wcs.utils.proj_plane_pixel_area(wcs_LW) ** 0.5 * 3600.
8 spectra files in total.
In [ ]:
 
In [92]:
# coords_bound = SkyCoord(['11:06:00.7328 -77:21:29.267', '11:06:51.4474 -77:22:46.359'], unit = (u.h, u.deg))
# opt_bound = np.polyfit(coords_bound.ra.value, coords_bound.dec.value, deg = 1)

for i, tmp_spec_path in enumerate(list_spec[1:2]):
    tmp_spec_fits = fits.open(tmp_spec_path)
    '''source information '''
    tmp_id = tmp_spec_fits[0].header['ID']         # tmp_spec_fits[0].header['NUMBER']
    tmp_filter = tmp_spec_fits[0].header['FILTER'] # tmp_spec_fits[0].header['select'] #
    tmp_mag = tmp_spec_fits[0].header['ABMAG']     # tmp_spec_fits[0].header['MAG_AUTO']
    # if tmp_mag < 23.5: continue
    tmp_flux_mJy = 10**(-0.4 * tmp_spec_fits[0].header['ABMAG']) * 3631e3 # 10**(-0.4 * tmp_spec_fits[0].header['MAG_AUTO']) * 3631e3 
    tmp_N_R = tmp_spec_fits[0].header['N_R']
    tmp_N_C = tmp_spec_fits[0].header['N_C']
    tmp_RA, tmp_DEC = tmp_spec_fits[0].header['RA0'], tmp_spec_fits[0].header['DEC0']
    tmp_x0, tmp_y0 = tmp_spec_fits[0].header['xcentroid'] , tmp_spec_fits[0].header['ycentroid'] 
    # tmp_x0, tmp_y0 = tmp_spec_fits[0].header['X_IMAGE'] - 1, tmp_spec_fits[0].header['Y_IMAGE'] - 1
    
    # if np.polyval(opt_bound, tmp_RA) > tmp_DEC: tmp_module = 'A'
    # else: tmp_module = 'B'
    
    '''2d spec'''
    tmp_spec_2d = tmp_spec_fits[1].data
    tmp_wht_2d = tmp_spec_fits[2].data
    wave_sample_c = tmp_spec_fits[1].header['WAVE_1'] + np.arange(tmp_spec_fits[1].header['NAXIS1']) * tmp_spec_fits[1].header['D_WAVE']
    
    ## if not enough wavelength coverage, then skip
    if np.sum(np.isnan(np.nansum(tmp_spec_2d, axis = 0)) == False) < 200: continue
        
    ## get y center:
    tmp_spec_ydir = np.nansum(tmp_spec_2d * tmp_wht_2d, axis = 1) / np.nansum(tmp_wht_2d, axis = 1)
    tmp_yc = np.argmax(tmp_spec_ydir[12:19]) + 12
    # tmp_yc = 15
    tmp_aper = 3
    
    ## model background with 1d fit at each x-pix position
    h_space_bkg = np.append(np.arange(0, tmp_yc-6), np.arange(tmp_yc+6, 31))
    tmp_specbkg_2d = tmp_spec_2d[h_space_bkg]
    _, tmp_med_specbkg_2d, tmp_rms_specbkg_2d = sigma_clipped_stats(tmp_specbkg_2d, axis = 0)
    tmp_specbkg_2d_3sigma = np.ones_like(tmp_specbkg_2d) * (tmp_med_specbkg_2d + tmp_rms_specbkg_2d * 2.5)
    tmp_wht_specbkg_2d = np.clip(np.float32(tmp_specbkg_2d < tmp_specbkg_2d_3sigma), 0.01, 1)
    tmp_specbkg_2d_model = np.vstack([np.polyval(np.polyfit(h_space_bkg, np.nan_to_num(tmp_specbkg_2d[:, x]), deg = 1, 
                                                            w = tmp_wht_specbkg_2d[:, x]),
                                                 np.arange(31)) for x in range(tmp_specbkg_2d.shape[1])])
    tmp_specbkg_2d_model = tmp_specbkg_2d_model.T
    tmp_spec_2d = tmp_spec_2d - np.nan_to_num(tmp_specbkg_2d_model)
    
    tmp_spec_1d = np.nansum(tmp_spec_2d[tmp_yc-tmp_aper:tmp_yc+tmp_aper+1], axis = 0) 
    tmp_unc_1d = np.nansum(tmp_wht_2d[tmp_yc-tmp_aper:tmp_yc+tmp_aper]**-1, axis = 0)**0.5
    tmp_unc_1d = ((tmp_rms_specbkg_2d / 4.)**2 + tmp_unc_1d**2)**0.5


    
    '''
    Plot Spectra:
    '''
    plt.close()
    fig, ax = plt.subplots(2, 1, figsize = (12, 5))
    
    ### ax[0]: 2D spectra
    ax[0].imshow(tmp_spec_fits[1].data[:], aspect = 5, vmax = np.clip(tmp_flux_mJy / 2., 0.01, 1e5), vmin = -0.01, 
                 cmap = plt.cm.gist_gray_r)
    ax[0].set(ylim = (0.5, 30.5), xticks = [], aspect = 6 * (np.diff(WRANGE) - 0.1) / 1.2,
              xlim = (0.05 / tmp_spec_fits[1].header['D_WAVE'], 
                      (WRANGE[1] - WRANGE[0] - 0.05) / tmp_spec_fits[1].header['D_WAVE'])
             )
    ax[0].set_yticks([15.5])
    ax[0].set_yticklabels([""])
    
    ### ax[1]: 1D spectra
    ax[1].plot(wave_sample_c, ndimage.gaussian_filter1d(tmp_spec_1d, 0.5), # tmp_spec_1d, #
               color = 'k', zorder = 100, lw = 1.5)
    ## uncomment below to add R/C-only spectra, if you have extracted them:
    # tmp_spec_R_path = tmp_spec_path.replace('_comb.fits', '_R.fits') 
    # tmp_spec_C_path = tmp_spec_path.replace('_comb.fits', '_C.fits')
    # if os.path.isfile(tmp_spec_R_path):
    #     tmp_r_fits = fits.open(tmp_spec_R_path)
    #     wave_sample_gR = tmp_r_fits[1].header['WAVE_1'] + np.arange(tmp_r_fits[1].header['NAXIS1']) * tmp_r_fits[1].header['D_WAVE']
    #     tmp_spec_R_1d = np.sum(tmp_r_fits[1].data[tmp_yc-3:tmp_yc+4], axis = 0) 
    #     tmp_specbkg_R_2d = np.vstack((tmp_r_fits[1].data[0:10], tmp_r_fits[1].data[20:]))
    #     tmp_bkg_R_1d = sigma_clipped_stats(tmp_specbkg_R_2d, axis = 0, sigma_upper = 2)[0]
    #     tmp_spec_R_1d = tmp_spec_R_1d - tmp_bkg_R_1d * 10
    #     ax[1].plot(wave_sample_gR, tmp_spec_R_1d , color = 'r', zorder = 50, lw = 0.5)
    # if os.path.isfile(tmp_spec_C_path):
    #     tmp_c_fits = fits.open(tmp_spec_C_path)
    #     wave_sample_gC = tmp_c_fits[1].header['WAVE_1'] + np.arange(0, tmp_c_fits[1].header['NAXIS1'], 1) * tmp_c_fits[1].header['D_WAVE']
    #     tmp_spec_C_1d = np.sum(tmp_c_fits[1].data[tmp_yc-3:tmp_yc+4], axis = 0) 
    #     tmp_specbkg_C_2d = np.vstack((tmp_c_fits[1].data[0:10], tmp_c_fits[1].data[20:]))
    #     tmp_bkg_C_1d = sigma_clipped_stats(tmp_specbkg_C_2d, axis = 0, sigma_upper = 2)[0]
    #     tmp_spec_C_1d = tmp_spec_C_1d - tmp_bkg_C_1d * 10
    #     ax[1].plot(wave_sample_gC, tmp_spec_C_1d , color = 'c', zorder = 49, lw = 0.8)
    tmp_max_counts = np.nanpercentile(tmp_spec_1d[(np.isnan(tmp_spec_1d) == False) & (tmp_spec_1d != 0)], 95) * 1.5
    if np.isnan(tmp_max_counts): continue
    
    ax[1].axhline(0, color = 'grey', ls = '--' )
    ax[1].set(xlim = (WRANGE[0] + 0.05, WRANGE[1] - 0.05), 
              xticks = np.arange(WRANGE[0] + 0.1, WRANGE[1] - 0.05 + 0.01, 0.1),
              xlabel = 'Observed Wavelength (µm)', ylabel = 'Flux Density [mJy]')
    ax[1].set_ylim(-0.03, np.clip(tmp_max_counts, 0.03, 1e8)) 
    
    corner_text(ax[0], loc = 2, s = 'ID%s' % (tmp_id), weight = 'semibold', color = 'r', fontsize = 18, edge = 5e-3)
    corner_text(ax[0], loc = 1, s = '%s=%.2f' % ('F430M', tmp_mag), color = 'r', fontsize = 14, edge = 5e-3)
    corner_text(ax[0], loc = 4, s = '(%.5f, %.5f)' % (tmp_RA, tmp_DEC), color = 'r', fontsize = 14, edge = 5e-3)
    corner_text(ax[1], loc = 1, s = '$f_\mathrm{%s}$=%.3f mJy' % ('F430M', tmp_flux_mJy), fontsize = 14, color = 'r', edge = 5e-3)
    corner_text(ax[1], loc = 2, s = 'N(R)=%2d' % tmp_N_R, fontsize = 12, color = 'r', edge = 5e-3)
    corner_text(ax[1], loc = 3, s = 'N(C)=%2d' % tmp_N_C, fontsize = 12, color = 'c', edge = 5e-3)
    fig.subplots_adjust(hspace = -0.20, wspace = 0, bottom = 0.13, top = 1.10, right = 0.98, left = 0.15)
    
    ### ax_im: Direct Image
    ax_im = fig.add_axes([0.02, 0.68, 0.30 / 12 * 5, 0.30])
    tmp_img = img_LW[int(tmp_y0)-16:int(tmp_y0)+17,int(tmp_x0)-16:int(tmp_x0)+17]
    ax_im.imshow(tmp_img, 
                 vmin = np.nanpercentile(tmp_img, 2.5), vmax = np.nanpercentile(tmp_img, 97.5),
                 cmap = plt.cm.gist_heat)
    ax_im.set(xlim = (tmp_x0%1 + 0.5, tmp_x0%1 + 30.5), ylim = (tmp_y0%1 + 0.5, tmp_y0%1 + 30.5))
    ax_im.set_xticks([])
    ax_im.set_yticks([])
    corner_text(ax_im, loc = 4, s = '(%.1f, %.1f)' % (tmp_x0, tmp_y0), color = 'w', fontsize = 10)
    
    print(i, tmp_id, '%s=%.2f' % (tmp_filter, tmp_mag) )
    tmp_spec_fits.close()
    '''
    Save spectral plot (remember to change the path)
    '''
    plt.savefig('/home/u24/fengwusun/IceAge/extracted_spectra/spec_1d_%s_%s_ID%s_mag%.2f.pdf' % (tmp_target, tmp_filter,
                                                                                                   tmp_id, tmp_mag), 
                dpi = 150)
    
    '''
    Save 1D spectra (remember to change the path)
    '''
    tb_1dspec = Table(names = ['wavelength_um', 'flux_mJy', 'fluxerr_mJy'], 
                      data = [wave_sample_c, tmp_spec_1d, tmp_unc_1d])
    tb_1dspec['wavelength_um'].info.format = '.4f'
    tb_1dspec['flux_mJy'].info.format = '.5f'
    tb_1dspec['fluxerr_mJy'].info.format = '.5f'
    ascii.write(tb_1dspec, tmp_spec_path.replace('/spec_2d_', '/spec_1d_').replace('.fits', '.dat'), 
                format = 'commented_header', overwrite = True)
2022-10-30 23:31:33,646 - stpipe - WARNING - /tmp/ipykernel_16422/670445884.py:51: RuntimeWarning: divide by zero encountered in reciprocal
  tmp_unc_1d = np.nansum(tmp_wht_2d[tmp_yc-tmp_aper:tmp_yc+tmp_aper]**-1, axis = 0)**0.5

0 1619 F444W=17.74

Optional: Check individual source & multiple frames

You can overlie spectral tracing / dispersion models on different 2D grism frames & individual sources to check the accuracy the spectral tracing models, and apply shifts for the POM catalogs if necessary.

In [102]:
'''
Select frames for F444W grism images, plot tracing/dispersion models for top-20 bright sources:
'''
#  glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[5-6]00*_nrc[a-b]long_rate_lv1.5.fits')
for k, tmp_rate_path in enumerate(all_rate_list):
    if k != 1: continue
    ### POM-applied catalog:
    tmp_tb_sl = ascii.read('/home/u24/fengwusun/IceAge/POM_catalog/%s_%s_dirimg_sources.list' % (os.path.basename(tmp_rate_path)[:-20], tmp_filter))
    ### LW grism image:
    tmp_rate_fits = fits.open(tmp_rate_path)
    tmp_rate_wcs_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
    tmp_grism_hd = tmp_rate_fits[0].header
    tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
    print('[%-3d]' % k, tmp_rate_wcs_base)
    print('/filter:', tmp_filter, '/module:', tmp_module, '/grism:', tmp_pupil)

    ### Source catalog - selected the brightest sources
    tmp_tb_sl.sort('MAG_AUTO')
    tmp_tb_sl = tmp_tb_sl[:20]
    

    ### Wave space for plotting
    if tmp_filter == 'F322W2': wave_space = np.arange(2.4, 4.11, 0.2)
    elif tmp_filter == 'F444W': wave_space = np.arange(3.85, 5.11, 0.2)

    ### Spectral tracing & Dispersion parameters:
    if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
    else: disp_filter = tmp_filter
    tb_order23_fit = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, tmp_module, tmp_pupil))
    fit_opt, fit_err = tb_order23_fit['col0'].data, tb_order23_fit['col1'].data
    tb_fit_displ = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % (tmp_module, tmp_pupil))
    w_opt, w_err = tb_fit_displ['col0'], tb_fit_displ['col1']
    
    '''Plot 2D grism spectra'''
    plt.close()
    fig, ax  = plt.subplots(1, 1, figsize = (12, 12))
    plt.imshow(tmp_rate_fits[1].data, vmin = -0.03, vmax = np.percentile(tmp_rate_fits[1].data, 98.), 
               origin = 'lower', aspect = 1, cmap = plt.cm.viridis)
    ### show stars
    for i, item_sl in enumerate(tmp_tb_sl):
        if item_sl['Index'] == 1: kwargs_target = dict(marker = '*', color = 'pink', mec = 'w', mew = 2, ms = 20)
        else: 
            kwargs_target = dict(marker = '*', color = 'r', mec = 'pink', mew = 2, ms = 20)
        plt.plot(item_sl['pixel_x'], item_sl['pixel_y'], ls = 'none', **kwargs_target)

        dx_space = func_fit_wave(np.vstack((item_sl['pixel_x'] * np.ones_like(wave_space), 
                                            item_sl['pixel_y'] * np.ones_like(wave_space), wave_space)), *w_opt)
        dy_space = fit_disp_order23(np.vstack((item_sl['pixel_x'] * np.ones_like(dx_space), 
                                               item_sl['pixel_y'] * np.ones_like(dx_space), dx_space)), *fit_opt)
        if tmp_pupil == 'R':
            plt.plot(item_sl['pixel_x'] + dx_space, item_sl['pixel_y'] + dy_space,
                     color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
            for j, tmp_dx in enumerate(dx_space):
                if ((item_sl['pixel_x'] + dx_space[j] > 2042) | (item_sl['pixel_x'] + dx_space[j] < 5)|
                    (item_sl['pixel_y'] + dy_space[j] > 2042) | (item_sl['pixel_y'] + dy_space[j] < 5)): continue
                plt.text(item_sl['pixel_x'] + dx_space[j], item_sl['pixel_y'] + dy_space[j] - 3 ,
                         s = '%.1f' % wave_space[j], color = 'r', ha = 'center', va = 'top', fontsize = 10)

        elif tmp_pupil == 'C':
            plt.plot(item_sl['pixel_x'] + dy_space, item_sl['pixel_y'] + dx_space,
                     color = 'r', lw = 0.8, marker = '_', ms = 10, ls = ':')
            for j, tmp_dx in enumerate(dx_space):
                if ((item_sl['pixel_y'] + dx_space[j] > 2042) | (item_sl['pixel_y'] + dx_space[j] < 5)|
                    (item_sl['pixel_x'] + dy_space[j] > 2042) | (item_sl['pixel_x'] + dy_space[j] < 5)): continue
                plt.text(item_sl['pixel_x'] + dy_space[j], item_sl['pixel_y'] + dx_space[j] - 3 ,
                         s = '%.2f' % wave_space[j], color = 'r', ha = 'right', va = 'center', fontsize = 10)

    ax.set(xlim = (0, 2048), ylim = (0, 2048))

    corner_text(ax, loc = 3, s = '%s mod%s (%s)' % (tmp_filter, tmp_module, tmp_pupil), 
                color = plt.cm.Reds(0.05), weight = 'bold', fontsize = 28, edge = 1e-2)
    
    plt.tight_layout()
    ### Uncomment below if you want to save the plot; remember to change the path.
    # fig.savefig('/home/u24/fengwusun/IceAge/tracing/%s_%s_%s_%s.pdf' % (tmp_rate_wcs_base, tmp_filter, tmp_module, tmp_pupil),
    #             dpi = 100)
[1  ] jw01309025001_02101_00001_nrcb
/filter: F444W /module: B /grism: R
In [ ]:
 
In [104]:
'''
Select one frame of F444W grism images, plot tracing/dispersion models for one "random" bright source:
'''

k = 1
tmp_rate_path = all_rate_list[k]

### POM-applied catalog
tmp_tb_sl = ascii.read('/home/u24/fengwusun/IceAge/POM_catalog/%s_%s_dirimg_sources.list' % (os.path.basename(tmp_rate_path)[:-20], tmp_filter))

### 2D grism image
tmp_rate_fits = fits.open(tmp_rate_path)
tmp_rate_wcs_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
tmp_grism_hd = tmp_rate_fits[0].header
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
print('[%-3d]' % k, tmp_rate_wcs_base)
print('/filter:', tmp_filter, '/module:', tmp_module, '/grism:', tmp_pupil)

### Source catalog - sorted by brightness, selected ones with relatively complete spectral coverage
tmp_tb_sl.sort('MAG_AUTO')
tmp_tb_sl = tmp_tb_sl[tmp_tb_sl['MAG_AUTO'] < 20.5]
if tmp_filter == 'F322W2':
    if tmp_pupil == 'C':
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_y'] - 1600) < 500)][:]
    elif (tmp_module, tmp_pupil) == ('A', 'R'):
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 1700) < 500)][:]
    elif (tmp_module, tmp_pupil) == ('B', 'R'):
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 400) < 500)][:]
elif tmp_filter == 'F444W':
    if tmp_pupil == 'C':
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_y'] - 700) < 600)][:]
    elif (tmp_module, tmp_pupil) == ('A', 'R'):
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 700) < 600)][:]
    elif (tmp_module, tmp_pupil) == ('B', 'R'):
        tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 1300) < 600)][:]

tmp_tb_sl = tmp_tb_sl[:50]
### random source, or give an ID
i_random = np.random.randint(len(tmp_tb_sl))
# i_random = 62
tmp_tb_sl = tmp_tb_sl[i_random:i_random+1]


### Wavelength space
if tmp_filter == 'F322W2': wave_space = np.arange(2.4, 4.11, 0.2)
elif tmp_filter == 'F444W': wave_space = np.arange(3.85, 5.11, 0.2)

### Spectral tracing and dispersion parameters:
if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
else: disp_filter = tmp_filter
tb_order23_fit = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, tmp_module, tmp_pupil))
fit_opt, fit_err = tb_order23_fit['col0'].data, tb_order23_fit['col1'].data
tb_fit_displ = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % (tmp_module, tmp_pupil))
w_opt, w_err = tb_fit_displ['col0'], tb_fit_displ['col1']

'''Plot 2D spectra cutout from the 2d grism image'''
plt.close()
fig, ax  = plt.subplots(1, 1, figsize = (12, 5))
if tmp_pupil == 'C':
    plt.imshow(tmp_rate_fits[1].data.T, vmin = -0.03, vmax = np.percentile(tmp_rate_fits[1].data, 96.), 
               origin = 'lower', aspect = 1, cmap = plt.cm.viridis, interpolation = 'gaussian')
elif tmp_pupil == 'R':
    if (tmp_filter, tmp_module) == ('F322W2', 'B'): tmp_vmax = np.percentile(tmp_rate_fits[1].data, 80.)
    else: tmp_vmax = np.percentile(tmp_rate_fits[1].data, 96.)
    plt.imshow(tmp_rate_fits[1].data, vmin = -0.03, vmax = tmp_vmax, 
               origin = 'lower', aspect = 1, cmap = plt.cm.viridis, interpolation = 'gaussian')
### show stars
for i, item_sl in enumerate(tmp_tb_sl):
    if item_sl['Index'] == 1: kwargs_target = dict(marker = '*', color = 'pink', mec = 'w', mew = 2, ms = 20)
    else: 
        kwargs_target = dict(marker = '*', color = 'r', mec = 'pink', mew = 2, ms = 20)
    plt.plot(item_sl['pixel_x'], item_sl['pixel_y'], ls = 'none', **kwargs_target)

    dx_space = func_fit_wave(np.vstack((item_sl['pixel_x'] * np.ones_like(wave_space), 
                                        item_sl['pixel_y'] * np.ones_like(wave_space), wave_space)), *w_opt)
    dy_space = fit_disp_order23(np.vstack((item_sl['pixel_x'] * np.ones_like(dx_space), 
                                           item_sl['pixel_y'] * np.ones_like(dx_space), dx_space)), *fit_opt) 
    if tmp_pupil == 'R':
        plt.plot(item_sl['pixel_x'] + dx_space, item_sl['pixel_y'] + dy_space,
                 color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
        for j, tmp_dx in enumerate(dx_space):
            if ((item_sl['pixel_x'] + dx_space[j] > 2042) | (item_sl['pixel_x'] + dx_space[j] < 5)|
                (item_sl['pixel_y'] + dy_space[j] > 2042) | (item_sl['pixel_y'] + dy_space[j] < 5)): continue
            plt.text(item_sl['pixel_x'] + dx_space[j], item_sl['pixel_y'] + dy_space[j] - 3 ,
                     s = '%.1f' % wave_space[j], color = 'r', ha = 'center', va = 'top', fontsize = 10)

    elif tmp_pupil == 'C':
        plt.plot(item_sl['pixel_y'] + dx_space, item_sl['pixel_x'] + dy_space, # + 2.0, 
                 color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
        for j, tmp_dx in enumerate(dx_space):
            if ((item_sl['pixel_y'] + dx_space[j] > 2042) | (item_sl['pixel_y'] + dx_space[j] < 5)|
                (item_sl['pixel_x'] + dy_space[j] > 2042) | (item_sl['pixel_x'] + dy_space[j] < 5)): continue
            plt.text(item_sl['pixel_y'] + dx_space[j], item_sl['pixel_x'] + dy_space[j] - 3 ,
                     s = '%.2f' % wave_space[j], color = 'r', ha = 'right', va = 'center', fontsize = 10)
### xlim/ylim of the plot
if tmp_pupil == 'R':
    ax.set(xlim = (item_sl['pixel_x'] + np.min(dx_space) - 20, item_sl['pixel_x'] + np.max(dx_space) + 20 ), 
           ylim = (item_sl['pixel_y'] + np.median(dy_space) - 7.5, item_sl['pixel_y'] + np.median(dy_space) + 7.5 ), 
           aspect = 'auto')
elif tmp_pupil == 'C':
    ax.set(xlim = (item_sl['pixel_y'] + np.min(dx_space) - 20, item_sl['pixel_y'] + np.max(dx_space) + 20 ), 
           ylim = (item_sl['pixel_x'] + np.median(dy_space) - 7.5, item_sl['pixel_x'] + np.median(dy_space) + 7.5 ), 
           aspect = 'auto')

corner_text(ax, loc = 3, s = '%s mod%s (%s)' % (tmp_filter, tmp_module, tmp_pupil), 
            color = plt.cm.Reds(0.05), weight = 'bold', fontsize = 28, edge = 1e-2)
corner_text(ax, loc = 1, s = tmp_rate_wcs_base, color = 'w')
# plt.tight_layout()
[1  ] jw01309025001_02101_00001_nrcb
/filter: F444W /module: B /grism: R
Out[104]:
Text(1262.3794977135256, 1438.7053288047252, 'jw01309025001_02101_00001_nrcb')
In [105]:
tmp_tb_sl
Out[105]:
Table length=1
Indexradecpixel_xpixel_yMAG_AUTOMAGERR_AUTOFLAGS
int64float64float64float64float64float64float64int64
1296166.590421-77.39821148.5221441.03519.8570.0030
In [ ]:
 
In [ ]: